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Abstract 

We review the recent developments in the theory of normal, normal self-dual 
and general complex random matrices. The distribution and correlations of the 
eigenvalues at large scales are investigated in the large N limit. The 1/N expansion 
of the free energy is also discussed. Our basic tool is a specific Ward identity 
for correlation functions (the loop equation), which follows from invariance of the 
partition function under reparametrizations of the complex eigenvalues plane. The 
method for handling the loop equation requires the technique of boundary value 
problems in two dimensions and elements of the potential theory. As far as the 
physical significance of these models is concerned, we discuss, in some detail, the 
recently revealed applications to diffusion-controlled growth processes (e.g., to the 
Saffman- Taylor problem) and to the semiclassical behaviour of electronic blobs in 
the quantum Hall regime. 
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1 Introduction 



The subject matter of random matrix theory is a matrix whose entries are randomly 
distributed with some probability density. To put it another way, the theory deals with 
statistical ensembles of matrices. Given such an ensemble, one is typically interested in 
the distribution of eigenvalues and correlations between them as size of the matrices, N, 
tends to infinity. The distribution and correlation laws obtained in this way turn out to 
be common to objects and systems of very diverse nature. 

The area of applications of the random matrix theory in physics (and mathematics) 
is enormously vast. It ranges from energy levels statistics in nuclei to number theory, 
from quantum chaos to string theory. Most extensively employed and best-understood 
are ensembles of hermitian or unitary matrices, with eigenvalues being confined either to 
the real line or to the unit circle. Their applications to the level statistics in nuclei go 
back to Wigner's works of early 50-s. For different aspects of random matrix theory, its 
applications and related topics see e.g. [l]-[4]. 

In these lectures we consider more general classes of random matrices, with no a 
priori restrictions to their eigenvalues being imposed. The eigenvalues can be arbitrary 
complex numbers. Such models are as yet less well understood but they are equally 
interesting and meaningful. As we shall see, they may exhibit even richer mathematical 
structures than their Hermitian counterparts. Their physical applications are also many 
and varied. (A list of the relevant physical problems and corresponding references can 
be found in, e.g., [5].) The present lectures are based on our recent works [6]- [10] where 
new applications to diffusion limited growth processes, complex analysis and quantum 
Hall effect were found. 

The progenitor of ensembles of matrices with general complex eigenvalues is the statis- 
tical model of complex matrices with the Gaussian weight. It was introduced by Ginibre 
[11] in 1965. The partition function of this model is 

Z N = J [£>$] exp ^- ytr &<5>^J 

Here [D<&] = Yl- d(lZe$ij)d(Im&ij) is the standard volume element in the space of 
N x N matrices with complex entries and t is a (real positive) parameter. Along 
with the Ginibre ensemble and its generalizations we also consider ensembles of normal 
matrices, i.e., such that $ commutes with its hermitian conjugate <&t, and normal self- 
dual matrices (the definition follows below in Section 2). 

Since one is primarily interested in statistics of eigenvalues, it is natural to express 
the probability density in terms of complex eigenvalues Zj = X{ + iyj of the matrix $. It 
appears that the volume element can be represented as 

[D^ocJllzi-Zj^llcPzi 

i<j i 

where (5 — 1 for complex and normal matrices and j3 = 2 for normal self-dual matrices. If 
the statistical weight depends on the eigenvalues only, as it is usually assumed, the other 
parameters of the matrix (often referred to as "angular variables" ) are irrelevant and can 
be integrated out giving an overall normalization factor. In this case the original matrix 
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problem reduces to statistical mechanics of iV particles with complex coordinates Zj in the 
plane. We thus see that even if the matrix entries $y are statistically independent, like 
in the Ginibre ensemble, the eigenvalues are correlated in a nontrivial way. Specifically, 
the factor Yli<j \ z i ~ z j\ ■> being equal to the exponentiated Coulomb energy in two 
dimensions, means an effective "repelling" of eigenvalues. This remark leads to the 
Dyson logarithmic gas interpretation [12], which treats the matrix ensemble as a two- 
dimensional "plasma" of eigenvalues in an external electric field. 

At (5 — 1, there is another important interpretation. Namely, the factor Yli<j ( z i ~ z j) 
can be thought of as coming from the Slater determinant of one-particle fermionic states. 
The averaging over matrices then turns into averaging over the ground state of a system 
of TV non-interacting fermions. In the case of the Ginibre ensemble it is the system of N 
electrons in a uniform magnetic field at the lowest Landau level. In the case of a spin-i 
electron in a non-uniform magnetic field all energy levels split, with the only exception 
of the lowest one, which remains highly degenerate. We shall see that the normal and 
complex matrix ensembles with a non-Gaussian statistical weight are equivalent to N 
polarized electrons in a non-uniform magnetic field confined to the lowest energy level. 
If the degeneracy of the level equals N, i.e., if the level is completely filled, the system 
of N electrons behaves as an incompressible quantum Hall droplet [13]. 

When N becomes large some new features emerge, which require a different language 
for their adequate description, in much the same way as classical thermodynamics results 
from statistical mechanics. As N — > oo, the eigenvalues densely fill a domain in the 
complex plane with the mean density outside it being exponentially small in N. Around 
the edge of this domain the density steeply drops down. The width of the transition 
region tends to zero as N — > oo, so that the density profile in the direction normal to 
the edge looks like a step function. This fact allows one to introduce the support of 
eigenvalues to be the region where the mean density of eigenvalues does not vanish as 
N — > oo. Typically, it is a bounded domain (or several disconnected domains) in the 
complex plane. Its shape is determined by the probability density. 

For the Ginibre ensemble, the support of eigenvalues is the disk of radius y/t with 
uniform density. It is the counterpart of the celebrated Wigner "semicircular law" . For 
matrix ensembles with non-Gaussian weights the supports are in general not circular 
and not connected. Throughout these lectures our attention is mostly restricted to the 
case when the support of eigenvalues is a connected domain. Even in this relatively 
simpler case, the shape of this domain depends on parameters of the statistical weight in 
a rather complicated way. As we shall see in Section 4, the problem to find the support of 
eigenvalues from a given statistical weight is equivalent to the inverse problem of potential 
theory in two dimensions. In most cases, solutions of the latter are not available in an 
explicit form. 

Nevertheless, the local dynamical law that governs the evolution of the support of 
eigenvalues under changes of parameters of the statistical weight (like t in the Ginibre 
ensemble) can be expressed in terms of the exterior Dirichlet boundary value problem. 
Namely, the edge of the support moves along gradient of a scalar harmonic field in 
its exterior, with the velocity being proportional to the absolute value of the gradient. 
Remarkably, this growth law is known to be common to a wide class of diffusion-limited 
growth processes of which the most popular example is viscous flow in the Hele-Shaw 
cell (see [14] for a review). The mentioned above equivalence between the normal matrix 
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ensemble and the quantum Hall droplet suggests that the semiclassical behaviour of 
electronic droplets in a non-uniform magnetic field follows the same laws as the Hele- 
Shaw flows do. 

This fact allows one to treat the model of normal or complex random matrices as a 
growth problem. The advantage of this viewpoint is two-fold. First, the hydrodynamic 
interpretation makes some of the large N matrix model results more illuminating and 
intuitively accessible. Second and most important, the matrix model perspective may 
help to suggest new approaches to the long-standing growth problems. In this respect, 
of special interest is the identification of finite time singularities in some exact solutions 
to the Hele-Shaw flows with critical points of the normal and complex matrix models. 

At last, a few words about the organization of the lectures. The material that follows 
can be divided into three parts. The first one (Section 2) can be regarded as a continuation 
of the introduction. We define the main matrix ensembles to be considered and give their 
physical interpretations. The second part (Section 3) contains exact results valid at any 
finite N. We outline the integrable structure of the normal and complex models at f3 — 1 
(the Hirota relations for the partition function, the orthogonal polynomials technique and 
the Lax representation). In addition, we derive the exact relation between correlation 
functions of the eigenvalue densities (referred to as the loop equation) which holds for 
arbitrary values of (5. In the third part (Sections 4 and 5) we examine the large N limit 
of the models of random matrices with complex eigenvalues and discuss the applications 
to the growth processes and to the semiclassical electronic droplets in magnetic field. 
The Appendices contain technical details of some proofs and calculations. 

2 Some ensembles of random matrices with complex 
eigenvalues 

We consider square random matrices $ of size N with complex entries $jj subject to 
certain constraints depending on the particular ensemble. Some ensembles of random 
matrices are listed in the following table: 



Ensemble 


Notation 


Condition 


Dimension 


Hermitian 


"K 


$t = $ 


N 2 


Unitary 


IX 


$t$ = i 


N 2 


Normal 


N 


[$t,$] = 


N 2 + N 


Normal 
self-dual 




= 
$ self-dual 


±N 2 + N 


Complex 


e 


none 


2N 2 
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The first two matrix ensembles, % and It, are the most popular ones. They are given 
here just for comparison. Eigenvalues of matrices from IK and IX are confined to the real 
axis and to the unit circle respectively. The last three ensembles (which are the main 
subject of these lectures) do not imply any a priori restrictions on eigenvalues of the 
matrices. Normal matrices are defined by the constraint that they commute with their 
adjoint. The ensemble N° is defined for N even only. The meaning of the condition "$ 
is self-dual" is explained below in this section. By dimension of the ensemble we mean 
the real dimension of the matrix variety. 

Throughout this paper we consider the probability densities of the form P($) oc 
e trW ^\ where the function W(<&) (often called the potential of the matrix model) is a 
matrix-valued function of $ and $^ such that (W(<&))t = W(<&). This form is similar to 
the one usually employed in Hermitian and unitary ensembles. The partition function is 
defined as the integral over matrices from one or another ensemble: 

Z N = J [D$]e trf « (2.1) 

Summing over N with a suitable weight, one may also define the grand canonical ensem- 
bles corresponding to (2.1) but we do not pursue this possibility here. 

We need to specify the integration measure [-D<3>] and the potential W(<&). Given the 
measure and the potential, one is usually interested in the distribution and correlations 
of the eigenvalues. In general, they can be distributed on the real line for "K, on the unit 
circle for It and on the whole complex plane for N, N° and C. 



2.1 Integration measures 

The integration measure has the most simple form for the ensemble of general complex 
matrices: 

N 

[£>$] = Y[ d{'Re<5> ij )d{Im<5> ij ) 

This measure is additively invariant and multiplicatively covariant, i.e. for any fixed 
(nondegenerate) matrix A e C we have the properties + A)] = [_D<J»] and [_D(<M.)] = 
[£)(A$)] = | det A| 2Ar [£>$]. The first one is obvious, to prove the second one is an easy 
exercise. It is clear that the measure is invariant under transformations of the form 
$ — > W&U with a unitary matrix U ("rotations" in the matrix space). 

The measure for N is induced by the standard flat metric in C, 

\\5<f>\\ 2 = tT(6$6&) = J2\^ij\ 2 

ij 

via the embedding N C C. Here N is regarded as a hypersurface in C defined by the 
quadratic relations $$t = The measure for the ensemble IN" is defined in a similar 
way. 

As usual in matrix models, we would like to integrate out the "angular" variables and 
to express the integration measure through eigenvalues of the matrices. 
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The measure for 3\f through eigenvalues [1, 15]. We derive the explicit represen- 
tation of the measure in terms of eigenvalues in three steps: 

1. Introduce coordinates in X C C. 

2. Compute the inherited metric on 3Sf in these coordinates: ||5$|| 2 = Qa^d^d^. 

3. Compute the volume element = sf\ det g af 3\ Y[ a d£ a . 

Step 1: Coordinates in N. For any matrix $, the matrices Hi = |($ + H 2 = 
^(<3> - $*) are Hermitian. The condition [$, $ f ] = is equivalent to [H 1: H 2 ] = 0. Thus 
Hi t 2 can be simultaneously diagonalized by a unitary matrix U : 

H 1 = UXW, X = dmg{x u ...,x N } 
H 2 = UYW, Y = diag{ yi ,...,y N ] 

Introduce the diagonal matrices Z = X + iY, Z = X — iY with diagonal elements 
Zj = Xj + iyj and Zj = Xj —iyj respectively. Note that Zj are eigenvalues of $. Therefore, 
any $ G N can be represented as 

$ = UZU* 

where U is a unitary matrix and Z is the diagonal matrix with eigenvalues of <3> on the 
diagonal. In fact normal matrices can be equivalently defined by the property of being 
the most general matrices that can be diagonalized by a unitary transformation. The 
matrix U is defined up to multiplication by a diagonal unitary matrix from the right: 
U — > U C/diag- The dimension of N is thus 

dim (N) = dim (IX) - dim (U diag ) + dim (e diag ) = N 2 - N + 2N = N 2 + N 

Let us make a remark that the naive counting of the number of constraints in the 
condition [<£>, <&t] — o leads to a wrong result for the dimension of N. On the first glance, 
this condition gives N 2 — 1 independent constraints. Indeed, set H = Then the 

conditions H tk = for I < k give N(N — 1) real constraints and the conditions H kk = 
give N — 1 real constraints (because tr if = identically), in total N 2 — 1 constraints. 
We thus observe that dim (Jsf) ^ dim(C) — (N 2 — l). Therefore, there are only N 2 — N 
independent constraints among the N 2 — 1 equations [<E>, <&t] = 0. This fact can be easily 
illustrated by the example of 2 x 2 matrices. 

Step 2: The induced metric. Since $ = U ZU\ the variation is 5<E> = U (5u ■ Z + 5Z + Z ■ 
5v))U\ where 8u^ = U5W = —5v). Therefore, 

||5$|| 2 = tr (5<Z>5&) = tr (5Z5Z) + 2 tr (5uZ5uZ-(5u) 2 ZZ) 

N N 
j=l j<k 

(Note that Sujj do not enter.) This is the square of the line element ||5$|| 2 = g a p8^ a 8^. 
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Step 3. The volume element. We see that the metric g a p is diagonal in the coordinates 
7Ze(5zj), Im(5zj), lZe(5iijk), Tm(5ujk) with 1 < j < k < N, so the determinant of the 
diagonal matrix g a/3 is easily calculated to be | det g a/3 \ = 2 N ~ N Y[j<k \ z i~ z k\ 4 - Therefore, 

N 

[D$] oc [DU]'\A N (z 1 ,...,z N )\ 2 l[d 2 z j (2.2) 

3=1 

where cPz = dxdy is the flat measure in the complex plane, [DU]' = [DU]/[DU < n s ,g\ is 
the invariant measure on lX/lXdi ag , and 

N 

A N ( Zl , ...,z N ) = U(zj-z k ) = det^z*- 1 ) (2.3) 
j>k 

is the Vandermonde determinant. 

Similarly to the ensemble "K of Hermitian matrices, the measure (2.2) contains the 
squared modulus of the Vandermonde determinant. The difference is that the eigenvalues 
are complex numbers. The statistical model of normal random matrices was studied in 
[15, 16]. 



Normal self-dual matrices. 

/ o 



Let T be the matrix 

1 \ 



1 
-1 



r 2 = -i 



V ••• / 

(all other entries are zero). A complex matrix $ is called self-dual if T$ T r = — $ (the 
superscript T means transposition). The size of a self-dual matrix is thus an even number. 
It can be shown that eigenvalues of self-dual matrices always come in pairs: the diagonal 
form of $ is 

Z = diag {zi, Zi, z 2 , z 2 , ■ ■ ■ , z^, z^} 

As is easy to verify, the condition that $ is self-dual is equivalent to the condition that 
the matrix r<£> is anti-symmetric. 

Normal self-dual matrices are parameterized as $ = U ZW with Z as above, where 
U is unitary and symplectic: WU = 1, U T TU = T. In other words, U belongs to the 
maximal compact subgroup in the complex group of symplectic matrices Sp(iV). The 
latter is known to have real dimension AN 2 + 2N ', with the dimension of the maximal 
compact subgroup being twice less. Therefore, similarly to the calculation for normal 
matrices, dim(N°) = 2N 2 + N-N + 2N = 2N 2 + 2N (the value given in the table above 
corresponds to N replaced by N/2). The integration measure appears to be 

N 

[£>$] oc [DU]'\A N (z 1 ,...,z N )\ 4 l[d 2 z j (2.4) 

i=i 

(for 2./V x 2N matrices). Note that the module of the Vandermonde determinant enters 
in the fourth degree. The statistical model of normal self-dual matrices was discussed in 
[17]. 
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The measure for C through eigenvalues. A complex matrix $ with eigenvalues 
Zi, . . . ,zn can be decomposed as 

$ = U(Z + R)U ] 

where Z = diag {z±, . . . , zjy} is diagonal, U is unitary, and R is strictly upper triangular, 
i.e., Rij — if % > j. These matrices are defined up to a "gauge transformation": 
U — > U E/diag, -R - ^ ^diag -R^diag- It is not so easy to see that the measure factorizes. 
This requires some work, of which the key step is a specific ordering of the independent 
variables. The final result is: 

[£>*] oc [DU]' (]Jd 2 R k A \^ N {z i )\ 2 \[d 2 z j ( 2 - 5 ) 
\k<i / 3=1 

The details can be found in the Mehta book [1]. 



2.2 Potentials 

For the ensembles 3\f, 74° the "angular variables" (parameters of the unitary matrix U) 
always decouple after taking the trace trVF(<E>) = £V W(zj), so the potential W can be 
a function of $, $^ of a general form W(<&) = a nm <& n (&) m . Two important particular 
cases arise if the potential is: 

• Axially symmetric, 1T( ( I ) ) = Wq($>&). In this case the iV-fold integral essentially 
reduces to ordinary ones, and so some basic results become available in a quite 
explicit form. 

• Harmonic on the background of $$ t , i.e., W{<&) = + V(<&) + V r ($ t ). In what 
follows, we call it quasiharmonic. Here V(z) is an analytic function of z in some 
domain containing the origin and V(z) = V(z). In terms of the eigenvalues, the 
quasiharmonic potential is 

W(z) = -\z\ 2 + V(z) + V{z) (2.6) 

This case is particularly important for applications. The normal matrix model 
with quasiharmonic potentials bears some formal similarities with the model of 
two coupled Hermitian matrices [18] and the matrix quantum mechanics in the 
singlet sector [19]. 

The partition function reduces to 

r N 

Z N = / |A^)rn^ (2 ^% (2.7) 

where f3 — 1 for K and (3 = 2 for N°. From now on this formula is taken as the 
definition of the partition function. Comparing to (2.1), we redefine W — > W/j3 and 
ignore a possible A^-dependent normalization factor. One may also consider this integral 
for arbitrary values of (5. 
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The choice of the potential for the ensemble C is more restricted. For a general 
potential, the matrix U in $ = U(Z + R)W still decouples but R does not. The problem 
becomes too complicated. An important particular case, when R nevertheless decouples 
is the quasiharmonic potential (2.6). Indeed, tr ($$"'") = tr (ZZ) + tr (RR^), tr ($ n ) = 
tr(Z + R) n = trZ n , and so 

jf[D$] e tr^W =Cn J \A( Zi )\ 2 l[e w ^d 2 z k (2.8) 

where Cat is an A^-dependent normalization factor proportional to the gaussian integral 
f{DR]e- tT ( RRt \ 

As an example, let us consider the quadratic potential: 

W(z) = -a\z\ 2 + 2TZe{t l z + t 2 z 2 ) , a > 

The ensemble 6 (f3 — 1) with this potential is known as the Ginibre-Girko ensemble 
[11, 20]. In this case the partition function (2.7) can be calculated exactly [21]: 

Z N = Z y N \o - 4|t 2 | ) / exp \N a 2_ A \ t2 \2 ) 



where 



N 



k=l 

To the best of our knowledge, there are no exact results for 2D integrals of this type with 
\A(zi)\ 213 for other values of (3, even for the pure Gaussian weight. 

Coming back to the general case, we note that some integrals considered above may 
diverge. In the most important case of quasiharmonic potential, the integral 



e ~\z\*+V(z)+V(z) d 2 z 

converges for potentials V(z) = az 2 + (3z + ^/Xjlog(z — a^) with |a| < | and fij > —1 
but it always diverges if V(z) is a polynomial of degree > 3. As usual in matrix models, 
really interesting science begins when integrals diverge! Let us say a few words about 
how one should understand divergent integrals. The conventional viewpoint is to treat 
all cubic and higher degree terms in the potential as a small perturbation. The integral 
is then regarded as a perturbative series for a theory which is believed to be well-defined 
on the nonperturbative level. The nonperturbative definition can be achieved either 
by introducing an ad hoc cutoff or via more sophisticated methods in the spirit of the 
Marinari-Parisi approach [22] (the "stochastic stabilization"). As far as the large A^ limit 
is concerned, the integral for the partition function can be defined through the expansion 
around the saddle point. This gives the l/A^-expansion 

log Z N ~ N~ h F {h) as N -> oo 
h 

Even if the integral for Z^ diverges, each term of the 1/N expansion is often well-defined. 
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2.3 Physical interpretations 



The ensembles of random matrices appear to be mathematically equivalent to some 
important model systems of statistical and quantum mechanics. They are: 

• The 2D Coulomb plasma (any (3) 

• Non-interacting fermions {(3 = 1) 

• Electrons in magnetic field ((3 = 1) 

The equivalence holds for any finite N. The first two interpretations are standard and 
well known. The third one can be regarded as a specification of the second one for models 
with eigenvalues distributed over the whole complex plane. Remarkably, in this very case 
the noninteracting fermions picture (which is rather formal for "K and II) acquires a very 
interesting physical content related to the quantum Hall effect. 

The Dyson gas picture. This interpretation, first suggested by Dyson [12] for the 
unitary, symplectic and orthogonal matrix ensembles, relies on rewriting |Aat(z;)| 2/3 as 

exp^X)^- 1°S \ z i~ z i\} ■ Clearly, the integral (2.7) looks then exactly as the partition 
function of the 2D Coulomb plasma (often called the Dyson gas) at "temperature" 1/(3, 
in the external electric field: 

Z N = J e -W*i,-,*N) Y[ d 2 z 3 (2.9) 
The eigenvalues play the role of the 2D Coulomb charges. The energy is 

E = - E lo § i* - z i i 2 - f 3 ' 1 E w (*s) ( 2 - 10 ) 

i<3 3 

The first sum is the Coulomb interaction energy, the second one is the energy due to 
the external field. For the ensembles % and IX the charges are confined to dimension 
1 (the real line or the unit circle) but interact as 2D Coulomb charges. So, the Dyson 
gas picture for the ensembles N, 3\f° and 6 looks even more natural. The Dyson gas 
interpretation becomes especially helpful in the large N limit, where it allows one to 
apply thermodynamical arguments. 

Non-interacting fermions ([3 = 1). Given any system of polynomials of the form 
P n (z) = z n + lower degrees, the Vandermonde determinant can be written as A N (zi) = 
det(^ fc_1 ) = det(Pfc_i(^j)). Let us rewrite the statistical weight as 

N 

\A N (z l )\ 2 l[e w ^ = \^(z 1 ,...,z N )\ 2 

3=1 

where ^(^i, . . . , z N ) = det^x^ (Pk-i(zj)e w ( z ^) is the (unnormalized) wave function of 
N non-interacting fermions (the Slater determinant), with one-particle wave functions 
being ipk(z) — Pk-i{z)e w ^ z \ The partition function is the normalization integral: 
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For the ensembles N and C (with quasiharmonic potential) the wave function \1> has a 
direct physical meaning as a wave function of 2D electrons in magnetic field. 



Electrons in the plane in magnetic field. Consider a charged particle with spin 
| (electron) moving in the plane in a strong (not necessarily uniform) magnetic field 
B = B(x,y) orthogonal to the plane. The Pauli Hamiltonian reads 



# = (UhV + A) 2 - ha 3 B) 
2m V / 



Here m is mass of the particle, 0-3 = diag(l, —1) is the Pauli matrix, B = d x A y — d y A x . 
In 2D, the complex notation is convenient: A = A x — iA y , A = A x + iA y , so that 
B = i(&A-dA). 

For a uniform magnetic field, B = B , one can choose the gauge A — ^-z. Solving 
the Schrodinger equation Hip = Eip, one gets Landau levels: 

HB ( 1 \ 1 

E n = [n+--s , n = 0,1,2,..., s = ±- 

m \ 2 J 2 

The gap between the levels is proportional to B /m. Each level is highly degenerate. 
The wave functions at the level E = are 

/ n f Bo, ,0 
i> n = z exp I —^\z\ 

We note in passing that the system of fermions in the magnetic field at the lowest energy 
level admits a collective field theory description which was discussed in the literature in 
different contexts (see e.g. [13, 23, 24]). 

Let us turn to the problem with a nonuniform magnetic field. Choose the gauge 
A = idW with a real- valued function W, then B = -2 ddW and div A = 8 A + dA = 0. 
Therefore, 

(tW + A) 2 = -Ah 2 dd + 2ih((dW)d + (dW)B) + \dW\ 2 
= (2hd + dW)(-2W + dW) - 2hddW 
The Hamiltonian can be represented as the 2x2 matrix 



H = 



H+ 
H 



where 2mH± = (iftV + A) 2 ± 2hddW. In general, the spectral problem for this Hamilto- 
nian does not admit an explicit solution. However, the level E = is very special. Note 
that H + factorizes: H + = {2hd + dW){— 2hd + dW), so exact wave functions at the level 
E — can be found by solving the first order equation 

H + i) = (2hd - dW)ip = 

The general solution is 

rl>(z) = P(z)exp^W(z)\ (2.11) 
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E = 

s = l/2 



s = -l/2 



Figure 1: The energy levels of a spin-| electron in a non-uniform magnetic field (schemat- 
ically). The lowest level E = remains highly degenerate. 

where P(z) is an arbitrary holomorphic polynomial. The zero energy level remains to 
be highly degenerate even in the nonuniform magnetic field (Fig. 1). This fact was first 
observed in [25]. 

To find degeneracy of the level, we solve the Poisson equation AW = —2B, 

W(z) = ~J log \z - (\B(()d 2 ( (2.12) 

and observe that W(z) tends to — - log|z| as \z\ — > oo, where 4> = J Bd 2 z is the total 
magnetic flux. If P(z) is of degree n, the asymptotics of ip for \z\ large is 

ij = P(z)e^ w{z) -> z n 

where O — 2vrfi, is the flux quantum. We require the wave functions to be normalizable, 
i.e, / \ijj\ 2 d 2 z < oo that means n < 0/0o — 1- Therefore, 

n max = [0/0o] - 1 

([...] is the integer part), and the degeneracy is equal to the number of flux quanta in 
the total flux: 

N = [0/0 o ] 

If the Coulomb forces can be ignored, the wave function of N electrons in the plane in the 
magnetic field at the lowest energy level is constructed as the N x N Slater determinant 
of the functions of the type (2.11) with polynomials of different degrees. 

The situation when the lowest energy level E = is completely filled, i.e., N = n max , 
is the (integer) quantum Hall (QH) regime. The notion of the QH droplet [13] implies 
that the electronic liquid is incompressible, i.e., all states at the lowest energy level are 
occupied. We come to the following conclusion: the QH droplet consisting of N electrons 
(in general, in a non-uniform magnetic field) is equivalent to the ensemble of normal 
N x N matrices. 
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3 Exact results at finite TV 



3.1 Correlation functions: general relations 

The main objects to be determined in random matrix models are correlation functions. 
In general, they are mean values of scalar- valued functions of matrices. The mean value 
of such a function of the matrix <£> is defined, with the help of the statistical weight, 
in the usual way: 

J [D$]A($)e tTW ^ 



(A) 



e trw(*) 



We shall consider functions that depend on eigenvalues only - for example, traces of 
matrices. Correspondingly, typical correlators which we are going to study are mean 
values of products of traces: (tr /(<&)), (tr /i(<3>) tr/ 2 (<&)) and so on. Clearly, they are 
represented as integrals over eigenvalues. For example, 



(tr /($)) = 



N f \A N (z l )\^f(z 1 )f[e w ^d 2 z, 
_ j=i 

r N 

/ lA^z^He^d'zj 

J 3=1 



Here, /(<&) = is any function of <3>t which is regarded as the function f(zi) = 

f(zi, Zj) of the complex argument Zi (and Zj) in the r.h.s. A particularly important 
example is the density function defined as 

p{z) = Y J ^z-z ] )=t,5{z-^) (3.1) 

3 

where S(z) is the two dimensional (^-function. As it immediately follows from the defini- 
tion, any correlator of traces is expressed through correlators of p: 

/n 
(p( Zl ) . . . p(z n )} h{ Zl ) . . . f n (z n ) J] d 2 z 3 (3.2) 

3=1 

Instead of correlations of density it is often convenient to consider correlations of the field 
<p(z) = -P lo S \ z - z i I' = -P lo S l det ( z ~ $ ) I' ( 3 - 3 ) 

3 

from which the correlations of density can be found by means of the relation 

4ir{3p(z) = -A<p(z) (3.4) 

Clearly, (p is the 2D Coulomb potential created by the eigenvalues (charges). As it directly 
follows from the definitions, 

13 



where (■ ■ -) N means the expectation value in the ensemble of N x N matrices. 

Handling with multipoint correlation functions, it is customary to pass to their con- 
nected parts. For example, in the case of 2-point functions, the connected correlation 
function is defined as 

(p{z 1 )p{z 2 )) c = (p{z 1 )p(z 2 )) - (p(z 1 )) (p(z 2 )) 

The connected multi-trace correlators are expressed through the connected density cor- 
relators by the same formula (3.2) with (p{zi) . . . p(z n )) c in the r.h.s. The connected part 
of the (n + l)-point density correlation function is given by the linear response of the 
n-point one to a small variation of the potential. More precisely, the following variational 
formulas hold true: 

5\ogZ N Hp(zi)} S 2 \ogZ N 

Connected multi-point correlators are higher variational derivatives of logZ^. These 
formulas follow from the fact that variation of the partition function over a general 
potential W inserts J2i$( z ~ z %) the integral. Let us stress that these formulas are 
exact for any finite N. 

For a later use, we mention the formula 

Mm )=ew(f:^((trm)) k )}j (3.6) 

which immediately follows from the expansion 

log (e tr A = log Z N (W + /) - log Z N (W) 



f 5\0gZ N 2 1 f 5 2 \0gZ N ,/>w/>A ,2 A >2 A / 

3.2 Integrable structure of the N and C ensembles ((3 = 1) 

The partition function (2.7) for (3 = 1, 



N 

Z N = / \A N (z t )\ 2 l[e w ^d 

J 3=1 



Z 3 



regarded as a function of iV and Taylor coefficients of the potential W , has remarkable 
properties, which we briefly review below. 



Determinant representation. The following simple but important determinant rep- 
resentation holds true: 

Z N = N\ det(Cij), l<hj<N (3.7) 
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where 

dj = J z i - x z'- x e w ^d 2 z (3.8) 

is the "matrix of moments" . The proof is almost a repetition of the corresponding proof 
for the hermitian model. Complexity of eigenvalues does not cause any difficulties. For 
completeness, the detailed proof is given in the Appendix. 



Hirota relations. Let us apply the determinant formula to (det(A — <&)): 

(det(A - $)) = J- / | Aat] 2 TT(A - Zj)e w ^#Zj = J- det [z i - l z^\\-z)e w d 2 z] 
Zn J . Z N 

3 

Comparing this with the determinant representation of Z N , we see that (det (A — $)) = 
X N Z~^ 1 det [0^ — A _1 Cj + ij]. Taking appropriate linear combinations of the lines, one can 
reduce this determinant to the determinant of a matrix which differs from C^ only in 
the last line. Similarly, in the determinant representation of (det(Ai — $) det(A2 — $)) 
only two last lines change. Now, some standard identities for determinants lead to the 
following relations: 

(Ai-A 2 ) (det(Ai-$) det(A 2 -$)) (det(A 3 -$)) + cyclic perm-s of (123) = (3.9) 

N 7 7 

<|det(A-$)| 2 ) 7V - |(det(A-$)) 7V | 2 = — ?±L^zl <|det(A-$)| Vi (3.10) 
Write 

where 

W + [A] = W(z) + log (l - |) 

is the potential modified by the (complex and multi- valued) term log (l — f ). (Since this 
term is always under the exp-function, there is no ambiguity in the choice of its branch.) 
In this notation, the above identities for determinants acquire the form of the Hirota 
bilinear equations [26]: 

(A1-A2) Z N (W+[X 1 ] + [X 2 \) Z N (W+[X 3 \) + cyclic perm-s of 1, 2, 3 = (3.11) 

Z N (W)Z N (W+[X] + [X}) - Z N (W+[X])Z N (W+[X}) 



+ 

Let us parameterize the potential as 



JL-\x\-^z N+1 (w)z N _ 1 (w+[x}+[x}) 



(3.12) 



k 

then these equations state that Z^/Nl, as a function of t k , t^, is the tau-function of the 
2D Toda lattice hierarchy. The transformation W — > W + [A] is equivalent to the change 
of variables t k — > t k — \X~ k which is known in the literature as the Miwa transformation. 
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Orthogonal polynomials and the kernel function. The orthogonal polynomials 
technique is useful not only for hermitian and unitary matrix ensembles but for the 
normal and complex ensembles as well. The orthogonal polynomials are introduced as 
mean values of the characteristic polynomials of the random matrices $: 

P B (A) = <det(A -<&))„ (3.13) 

Clearly, P n are polynomials in A of the form P n (X) = \ n + lower degrees. 

The main property of the polynomials introduced is their orthogonality in the complex 
plane: 

J P n (z)lUz)e w ^d 2 z = h n 5 mn (3.14) 
The square of the norm h n = ||-P n || 2 is connected with the partition function as 

h n — —-t -p^- , Z N = N\l[h n 
n+1 Z n 11 



n=0 



Again, the proof is completely parallel to the corresponding proof in the hermitian models. 
See Appendix for details. 

The functions 

V h n -i 

are orthonormal: 



I 



ip n (z)ip m {z)d z = 5 mn (3.15) 



These ?/> n 's are "one-particle wave functions" of electrons in the magnetic field. The 
iV-particle wave function is ^jv^i, • • • ,zn) ~ det^^fc)]. The joint probability to find 
"particles" at z±, . . . , zn is 

I^Mzi,...,**)! 2 = ^y|det[^-(^)]| 2 

Since |detM| 2 = det(MM f ), we can write: | det[^(^)]| 2 = det (j]n=i i>n(zj)i>n(zk)) ■ 
The expression under the determinant is called the kernel function: 

N 

K N (z, w) = J2Mz)Mw) (3.16) 

71=1 

The main properties of the kernel function are: 



• Hermiticity: Kn(z, w) — K N (w, z); 

• Normalization: J K N (z, z)d 2 z = N; 

• Projection property: f Kn(zi, z)Kn(z, z,2)d 2 z = Kn(zi, z 2 ). 
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(p(z 1 )p(z 2 )) N = 



+ K N (z 1 ,z 1 )5(z 1 -z 2 ) 



All density correlation functions can be expressed through the kernel function. For 
example: 

(p( z )) N = k n{z,z) 

K N (z 1 ,z 1 ) K N (z 1 ,z 2 ) 

K N (z 2 ,z 1 ) K N (z 2 ,z 2 ) 
The last term is a contact term. It does not contribute if z\ ^ z 2 . In general, one has: 
{p(zi) . . . p{z n )) N = det(K N (zi, Zj))i<i,j< n + contact terms 

where the contact terms vanish if all the points zi are different. 

Note that for the ensemble N with an axially-symmetric potential W(&) = W (&&) 
the orthogonal polynomials are simply P n (z) = z n and h n is given explicitly: 



The kernel function is 



/■oo 

lz \2n e Wo(\z\*) d 2 z = 2n / r 2n + l e W ^) dr 

Jo 



N-1 



K N (z,w) = e^ w ^ +w ^J2 



n=0 



{zw) n 

h n 



For example, for the Gaussian model with W(z) = —\z\ 2 the squared norms are h n = irn\ 
and the mean value of density is given by 



7T * — * 



N - 1 |_|2n 

2 \ ^ \Z\ 



n=0 



The Lax representation. The orthogonal polynomials obey the recurrence relation 
of the form zP n (z) = J2 k<n c nk P k (z). In terms of the -^-function it reads 

zip n (z) = r n ip n+1 (z) +^u k (n)ip n - k (z) 

k>0 

One can represent it as a "spectral problem" Lip = zip for the difference operator 

L = r n e 9 ^ + J2Mn)e- kd/dn (3.17) 

k>0 

(the Lax operator). Here e d l dn is the shift operator n — > n + 1 with the characteristic 
property e d / dn f{n) = f(n + l)e d / dn . 



Let 



W (z)^W^(z) + J2(hz k + t k z k ) 
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It can be shown that the dependence on the parameters tk, t k is given by the 2D Toda 
hierarchy 

d d 
—L = [A k , L] , -^L = [L, A k ] 
otk ot k 

where A k = (L k ) + + ±(L fc ) , A k = (Lt*)_ + and V = e~ d l dn r n + Y Jk > Q e kd l dn u k (n) 

is the conjugate Lax operator. Given an operator of the form O = Ylk b k e kd ^ dn , we use 

the standard definition (6)+ = T /k >o b k ekd/dn ^ (&)- = 'Ek<o b k ekd/dn and (O)o = b . 
The structure of the Toda hierarchy in models of random matrices was first revealed in 
[27], see also review [4]. 

In the case of quasiharmonic potential, it is convenient to modify the -^-functions: 

ipn -> Xn = 7= — e y(z) P n _i(2:) 

V hn-1 



The functions Xn obey the orthogonality condition: J Xn(z)Xm(z)e ' 2 ' d 2 z = 5 nm . Then 
we have two compatible linear problems: 

(Lx)n = ZXn , (L' t x)n = ^zXn 

The second equation can be proven by comparing the matrix elements of the both sides 
using integration by parts. 



3.3 The loop equation 

In matrix models, loop equations are exact relations which follow from the fact that the 
matrix integral defining the model does not depend on changes of integration variables. In 
our case, we may start directly from the integral over eigenvalues (2.7), thereby extending 
the result to any value of j3. 

Clearly, the integral (2.7) remains the same if we change the integration variables 
Zi — > Zi. In other words, it is invariant under reparametrizations of the z-coordinate, 
which we write in the infinitesimal form as Zi — > Zi + e(zi), Zi — > Zi + e(zi). For the integral 
Z N = J 6 -PE(.zi,-,zn) (f z . £ given in (2.10), the reparametrization yields, in the 
first order: 



IP 



n^ 2 



„ (8E , . BE \ 
The invariance of the integral is then expressed by the identity 

E/|-(^.) e -' E )n^=° 



valid for any e. Introducing a suitable cutoff at infinity, if necessary, one sees that the 
2D integral over z^ can be transformed, by virtue of the Green theorem, into a contour 
integral around infinity and so it does vanish. 
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Let us take e(z;) = -3—, where 2; is a complex parameter. The singularity at the 
point z does not destroy the above identity since its contribution is proportional to the 
vanishing integral § dzi/(zi — z) over a small contour encircling z. Therefore, we have 
the equality 

E/[-flf^]--n^=o 

where d z .E = - V — /T 1 dW(zi) (see (2.10)). Using the identity 

rr 1 zi- zi 

^ — > 1 ^ — > 2 ^ — v 1 

we rewrite it in the form (T(zi, . . . , zn)) = 0, where 



l \ l / z 



This identity gives an exact relation between one- and two-point correlation functions. 
To see this, we rewrite it in terms of ip(z) = —f3 J2 i log \z — zi\ 2 using the rule £V f(zi) = 
f f(z)p(z)d 2 z. The result is the loop equation 



I / A = {m) (3,8) 



where 

T(z) = (d<p(z)y + (2-(3)dMz) (3-19) 
The correlator at coinciding points is understood as ((dip(z)) 2 ) = lim (dtp(z) dip(z')). 

z' — >z 

We have got an exact relation between one- and two-point correlation functions, valid 
for any finite N. For historical reasons, it is called the loop equation. One may read it 
as a Ward identity obeyed by correlation functions of the model. Being written in the 
form (3.18), (3.19), it resembles conformal Ward identities. Since correlation functions 
are variational derivatives of the free energy, the loop equation is an implicit functional 
relation for the free energy. However, it is not a closed relation. It can be made closed 
by some additional assumptions or approximations. A combination with 1/N expansion 
is particularly meaningful. 



4 Large TV limit 

Starting from this section, we study the large iV limit of the random matrix models 
introduced in Section 2. Our main tool is the loop equation. We shall see that in the 
large iV limit meaningful geometric and algebro-geometric structures emerge, as well as 
important applications in physics. 
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4.1 Preliminaries 



In order to be prepared for taking the large N ( "quasiclassical" ) limit, it is convenient 
to introduce the "Planck constant" h by the rescaling W(z) — > jW(z), so the integral 
for the partition function acquires the form 

Z N = J [D*]e& TW W oc j \A N \ 2 ?\{e^ w ^d 2 Z] (4.1) 

which is ready for an ^-expansion. 

Now we can specify what we mean by the large N limit. Namely, we are going to 
consider the integral (4.1) in the limit 

N -> oo, h->0, HN = t finite, 

where t is a (positive) parameter having the dimension of area, or, equivalently, the 
integral 

/• / v \ 

j2. 



Z N = J \A N (z t )\ 2 ?l[exp (jW(z s )j d' 



Z 3 



3 

as N — > oo. With this convention, the N — > oo and h — > limits mean the same. It is 
natural to expect that Zn — ° e' 1 2f °^\ where the rescaled free energy Fq is a smooth 
function of t (and parameters of the potential) finite as iV — > oo. 

From now on we change the normalization of the functions p and if (see (3.1), (3.3)) 
multiplying them by the h: 

p{z) =hJ2S(z- Zi ) , <p(z) = -h(3^og \z - Zl \ 2 (4.2) 

i i 

and use these definitions hereafter. (The idea is to make their mean values finite as 
N — > oo.) The relation (3.4) between these functions remains unchanged. The density 
is now normalized as follows: 

J p(z)d 2 z = t 

Note that in our units the h has dimension of [length] 2 , p(z) is dimensionless and the 
partition function defined by (4.1) has dimension [length] N (f 3N + 2 -l 3 )^ j. e-) the combination 

h-fr*+fr-W Z N = J \A N (zi/Vh) \ 2(3 f[ (e&M ^ (4.3) 

is dimensionless. 

In terms of the renormalized ip, the loop equation (3.18) has the form 

I / dw ^^ d 2 c + (w ,)f) + s (avw> = o (4.4) 

where 

e = (2-P)h (4.5) 

is small as h — > 0. Note that e is exactly zero for the ensemble N° ((3 = 2), and so the 
last term does not enter the loop equation in this case. It is convenient to treat e and H 
as independent small parameters. 



20 



4.2 Solution to the loop equation in the leading order 

It is instructive to think about the large N limit under consideration in terms of the Dyson 
gas picture. Then the limit we are interested in corresponds to a very low temperature 
of the gas, when fluctuations around equilibrium positions of the charges are negligible. 
The main contribution to the partition function then comes from a configuration, where 
the charges are "frozen" at their equilibrium positions. It is also important that the 
temperature tends to zero simultaneously with increasing the number of charges, so the 
plasma can be regarded as a continuous fluid at static equilibrium. In the noninteracting 
fermions picture, this limit has some features of the quasiclassical approximation. Math- 
ematically, all this means that the integral is evaluated by the saddle point method, with 
only the leading contribution being taken into account. As h — > 0, correlation functions 
take their "classical" values (ip{z)) = p d (z), and multipoint correlators factorize in the 
leading order: (dip (z) dip (z 1 )) = dip c i(z)dip c i(z'). Then the loop equation (4.4) becomes 
a closed relation for ip d : 

i_ J aw(OAMO d , ( + ( a ^ (2) ) 2 + EaVd(2) = (4 . 6) 

Note that we hold the last term which is apparently of the next order in h. The role of 
this term will be discussed below. 



The case (3 = 2 (the ensemble 3\f°). We begin with the case (3 = 2, when the last 
term in the r.h.s. of (4.6) vanishes exactly. Let us apply d to both sides of the equation. 
This yields: 

-dW{z)Aip d {z) + dip d {z)Aip d {z) = 
Since Ap d (z) oc p c i(z) (see (4.2)), we obtain 

p cl (z) [dip d (z)-dW(z)}=0 (4.7) 

This equation should be solved with the additional constraints f p d (z)d 2 z = t (normal- 
ization) and pd(z) > (positivity). The equation tells us that either dip d (z) = dW(z) 
or Pd(z) = 0. Applying d, we get Aip d (z) = AW(z). Since Ap = —An(3p, this gives the 
solution for p d : 

PM = -^P$- "in the bulk" (4.8) 
47rp 

Here, "in the bulk" just means "in the region where p c i > 0" . As we shall see below, this 
result holds true, up to some details, for other values of (3 as well, so (3 is kept in this 
formula and in some formulas below. 

The physical meaning of the equation dip c i(z) = dW(z) is clear. It is just the condition 
that the charges are in equilibrium (the saddle point for the integral). Indeed, the 
equation states that the total force experienced by a charge at any point z, where p d ^ 0, 
is zero. The interaction with the other charges, dipd(z), is compensated by the force 
dW(z) due to the external field. 
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Figure 2: The support of eigenvalues. 



Support of eigenvalues. Let us assume that 




1 



AW(z) > 



(4.9) 



An 



For quasiharmonic potentials, cr(z) = 1/tt. If, according to (4.8), p c i = a j 1 (3 everywhere, 
the normalization condition for p c i can not be satisfied! So we conclude that p c i = a/ (3 in 
a compact bounded domain (or domains) only, and outside this domain one should switch 
to the other solution of (4.7), p c \ = 0. The domain D where p c i > is called support 
of eigenvalues (Fig. 2). In general, it may consist of several disconnected components. 
The complement to the support of eigenvalues, D c = C \ D, is an unbounded domain in 
the complex plane. For quasiharmonic potentials, the result is especially simple: p c i is 
constant in D and in D c . 

In terms of the mean value of the function ip(z) in the leading order, the above result 
reads 



As it follows from the theory of potential in two dimensions, this function is continuous 
across the boundary of D together with its first derivatives. However, the second order 
derivatives of this function have a jump across the boundary. 

To find the shape of D is a much more challenging problem. It appears to be equivalent 
to the inverse potential problem in two dimensions. The shape of D is determined by the 
condition dip c i(z) = dW(z) (imposed for all points z inside D) and by the normalization 
condition. Since d(f c i(z) = —(3 J pcl ^ i (see (4.2)), we write them in the form 



The integral over D in the first equation can be transformed to a contour integral by 
means of the Cauchy formula (see Appendix B). As a result, one obtains: 



This means that the domain D has the following property: the function dW(z) on its 
boundary is the boundary value of an analytic function in its complement D c . 




< 




for all z e D 




for all z eD. 



(4.10) 
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We continue our analysis for the quasiharmonic case, where 

dW{z) = -z + V\z) , AW(z) = AddW{z) = -4 
The normalization than means that the area of D is equal to flirt. Assume that: 

- V(z) = J2t k z k is regular in D (say a polynomial) 

- G D (it is always the case if —W has a local minimum at 0) 

- D is connected 

Then the first equation in (4.10) acquires the form 



/ = V\z) for z G D 

JdD C - * 



2m 



Expanding it near z = 0, we get: 

''-ikL^'-rd™ (4 ' u) 

We see that the "coupling constants" t k are harmonic moments of D c = C \ D and the 
area of D is nf3t. 

It is the subject of the inverse potential problem to reconstruct the domain from 
its area and harmonic moments. (In the case when the support of eigenvalues has sev- 
eral disconnected components, some additional conditions are required.) In general, the 
problem has many solutions. But it is known that locally, i.e., for a small enough change 
t — > t + 5t, t k — > t k + St k the solution is unique. 



The support of eigenvalues: a fine structure {(3 ^ 2). Let us take a closer look at 
the support of eigenvalues, taking into account the so far ignored term in (4.6). Applying 
d to the both sides of (4.6) yields, instead of (4.7), the equation 

dip cl {z) - dW{z) + £ -d\ogA Vcl {z) = 

with the extra term proportional to e (4.5). Acting by d once again, we obtain the 
equation for p d 

- ^A\ogp cl (z)+Pp cl (z) = a(z) (4.12) 

which looks like the Liouville equation in the "background" o~{z). The first term seems 
to be negligible as h — > 0. However, one should be careful since the small parameter 
stands in front of the term with the highest derivative. As a matter of fact, this term is 
negligible only if p d is not small! Indeed, in a region where p c \ ~ e~ N , the term eA log p c i 
is of order 1 and thus plays the dominant role. 

In fact the equation states that p c i never vanishes exactly but can be exponentially 
small as iV — > oo. In the bulk, where the charges are distributed with nonzero density 
at iV = oo, equation (4.12) systematically generates corrections to the value a(z)/j3. If 
a 7^ const, there are power-like corrections in e, as well as exponentially small ones. We 
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conclude that the effect of the extra term is negligible everywhere except the very vicinity 
of the edge of the support of eigenvalues. Therefore, the result for p d can still be written 
in the form 

p cl (z) =p- 1 a(z)e(z;D) (4.13) 

where Q(z; D) is the characteristic function of the domain D (which is 1 in D and in 
D c ). The role of the term eAlog p c \ is to make the edge smooth. Around the edge, the 
density rapidly (but smoothly) drops down to zero over distances of order y/e. So, the 
boundary has got a "fine structure" . 

All this is in agreement with the form of the first nonvanishing correction to the mean 
density found from the loop equation (4.4). Let us write (<p(z)) = (p c i{z) + <pn(z), where 
<ffi is of order h. The result for the ipn can be compactly written in terms of the function 

X (z) = log y/ncr(z) (4.14) 

and its harmonic continuation x H { z ) from the boundary of D to its exterior. Specifically, 
X H { Z ) is a harmonic function in D c (regular at oo) such that x H { z ) = x( z ) on the 
boundary. In other words, it is the solution of the (exterior) Dirichlet boundary value 
problem (see below). The loop equation yields 

( -e [ X ( z ) - X H (oo) + \] , zeD 
<p h (z) = { (4.15) 
( -e[x H (z)- X H (oo)], zeD c 

(Note the discontinuity of this function across the boundary.) The corresponding correc- 
tion to the mean density pn = (p) — p c \ is 

Pk(z) = ^ (e(z; D)Ax(z) - 5(z; dD)d n ( x (z) - X H ( Z )) ~ \#(z; dD)j (4.16) 

where S(z; dD) is the delta function with the support on the boundary and S'(z; dD) is 
its normal derivative (see Appendix B). Here and below, d n is the normal derivative at 
the boundary, with the normal vector being directed to the exterior of the domain D. 
The correction is so singular because the zeroth approximation (4.13) is singular by itself. 
The singular function p h is to be understood as being integrated with any smooth test 
function. The first term in (4.16) is a correction to the bulk density. The second one is 
a correction to the shape of the support of eigenvalues (it describes a small displacement 
of the edge). The third term signifies the presence of a double layer of charges around 
the boundary. This just means that the boundary is smoothed out. 

For the normal self-dual matrices (/3 = 2) the correction p% vanishes. Certainly, this 
does not mean that the boundary is sharp. It becomes smooth if higher corrections in h 
(caused by fluctuations of the particles) are taken into account. 

For the normal matrix model (j3 = 1) with quasiharmonic potential equation (4.12) 
reads ^ 

- -Alog p d (z) + 7tpd(z) = 1 
o 

The obvious solution is p c i = But it is not normalizable! One must look for another 
solution. The right solution differs from the constant by exponentially small terms in 
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<p> 




exponential decay 




Figure 3: The mean density profile for models with quasiharmonic potential in the large 
N limit. The correction to the pure step function is shown by the dotted line. The effect 
of the correction is to form a double layer of charges near the edge. 

the bulk but exhibits an abrupt drop across the boundary of the domain determined by 
the harmonic moments (4.11). So, up to exponentially small corrections, the solution is 
given by (4.13) with a — 1/ir and (3 = 1. It is instructive to note that the first two terms 
in (4.16) vanish at a — const but the third term does not, so the leading correction in 
the quasiharmonic case merely makes the edge smooth that results in a double layer of 
"charges" near the edge (see Fig. 3). 

From the support of eigenvalues to an algebraic curve. There is an interesting 
algebraic geometry behind the large N limit of matrix models. For simplicity, here we 
consider models with quasiharmonic potentials. 

In general, the boundary of the support of eigenvalues is a closed curve in the plane 
without self-intersections. The following important fact holds true. If V'(z) is a rational 
function, then this curve is a real section of a complex algebraic curve of finite genus. In 
fact, this curve encodes the 1/N expansion of the model. In the context of Hermitian 
2-matrix model such a curve was introduced and studied in [28, 29, 30]. 

To explain how the curve comes into play, we start from the equation d(p c i = dW, 
which can be written in the form z — V'(z) = G(z) for z G D, where 



Clearly, this function is analytic in D c . At the same time, V'(z) is analytic in D and all 
its singularities in D c are poles. Set 



Then S(z) = z on the boundary of the support of eigenvalues. So, S(z) is the analytic 
continuation of z away from the boundary. Assuming that poles of V are not too close to 
<9D, S(z) is well-defined at least in a piece of D c adjacent to the boundary. The complex 




S(z) 



V'(z) + G(z) 
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conjugation yields S(z) = z, so the function S(z) = S(z) must be inverse to the S(z): 

S(S(z)) = z 

("unitarity condition"). The function S(z) is called the Schwarz function [31]. 

Under our assumptions, the S(z) is an algebraic function, i.e., it obeys a polynomial 
equation R(z, S(z)) = of the form 

d+1 

R(z,S(z))= Y, a niz n {S{z)) 1 = 

n,l=l 

where = a n \ and d is the number of poles of V'(z) (counted with their multiplicities). 
Here is the sketch of proof. Consider the Riemann surface E = D c U <9D U (D c )* (the 
Schottky double of D c ). Here, (D c )* is another copy of D c , with the local coordinate z, 
attached to it along the boundary. On E, there exists an anti-holomorphic involution 
that interchanges the two copies of D c leaving the points of <9D fixed. The functions z 
and S(z) are analytically extendable to (D c )* as S(z) and z respectively. We have two 
meromorphic functions, each with d+1 poles, on a closed Riemann surface. Therefore, 
they are connected by a polynomial equation of degree d + 1 in each variable. Hermiticity 
of the coefficients follows from the unitarity condition. 

The polynomial equation R(z, z) = defines a complex curve r with anti-holomorphic 
involution (z, z) \— > (z, z). The real section is the set of points such that z — z. It is the 
boundary of the support of eigenvalues. 

It is important to note that for models with non-Gaussian weights (in particular, with 
polynomial potentials of degree greater than two) the curve has a number of singular 
points, although the Riemann surface E (the Schottky double) is smooth. Generically, 
these are double points, i.e., the points where the curve crosses itself. In our case, a 
double point is a point z^ E D c such that S(z^) = z^ but z^ does not belong to the 
boundary of D. Indeed, this condition means that two different points on £, connected 
by the antiholomorphic involution, are stuck together on the curve r, which means the 
self-intersection. The double points play the key role in deriving the nonperturbative 
(instanton) corrections to the large N matrix models results (see [32] for details). 

Finally, let us point out that a complex curve 

pin) 

can be associated to ensembles of 
finite matrices as well (at least for 13 = 1). For the model of two Hermitian matrices this 
was done in [33]. If the linear spectral problem for the L-operator is of finite order (see 
the end of Section 3.2), the curve can be defined as the "spectral curve" of the difference 
spectral problems (Lx) n — z Xn, (L^x)n = %Xn- Since the operators L and V< do not 
commute, the curve depends on n. (See [8] for details.) In contrast to the curve r the 
curve is in general a smooth curve. A properly performed n — > oo limit of this curve 
coincides with the complex curve r constructed from the support of eigenvalues. 

4.3 The free energy 

In this subsection we assume that the support of eigenvalues is connected. As is known, 
the free energy admits a 1/N expansion. We prefer to work with the equivalent h- 
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expansion, thus emphasizing its semiclassical nature. The first few terms of the H- 
expansion for the ensembles of normal matrices with a general potential are 



\ogZ N = c(N) + ^ + ^l + F 1 + 0(H) (4.17) 
ur a 

The explicit form of the c(N) is given below. In fact this term can be absorbed into a 
normalization. For example, one may normalize Zn dividing it by the partition function 
of the Gaussian model. Although in general the ft-expansion does not look like a topolog- 
ical one, it appears to be topological (i.e., only even powers of h enter) for the ensemble 
3\f° with arbitrary potential and for the ensemble N with quasiharmonic potential. 

The leading order. The partition function is given by (2.9), Z N = J e~ f5E Y\d 2 Zj, 
where 

-PE( Zl , . . . , z N ) = (3 log \zi - zj\ + 7T 1 W( Zj ) (4.18) 

i¥=3 3 

Writing the energy in terms of the density function, we have, in the leading order: 

-(3U 2 E\p\ = (3 j j p(z)p(() log \z - (\d 2 zd 2 ( + J W(z)p(z)d 2 z (4.19) 

We need to find the minimum of E[p] with the constraint j pd 2 z = t. This is achieved 
by variation of the functional E[p] + \(f pd 2 z — i) with the Lagrange multiplier A. The 
resulting equation is 



2(3 J log \z - (\p(()d 2 ( + W(z) + A = 



Upon taking the z-derivative, we see that the extremal p(z) is equal to the p c i(z), as 
expected, and the equation coincides with the previously derived one, dip c i(z) = dW(z), 
with 

Vd(z) = -(3 j log \z - (\ 2 p d (()d 2 ( = - J log \z - C| V(C)rf 2 C 

Assuming that W(0) = and D is connected, the Lagrange multiplier is fixed to be 
A = y?d(0)> an d so W(z) = <fci(z) — Vcz(O)- Plugging this into (4.19), we find the leading 
contribution to the free energy E /h 2 = max p (— (3E[p\) = —{3E\p c \\: 



F = ~ I I a(z) log 



1 1 

~z~( 



a(()d 2 zd 2 ( (4.20) 



which is basically the electrostatic energy of the domain D charged with the density a(z) 
with a point-like compensating charge at the origin. 

Since the t-derivative of the extremal value of the functional is equal to the Lagrange 
multiplier (with the sign minus), d t F = —A, we incidentally obtain the useful formula 

d t F = 2 log \z\cr(z)d 2 z (4.21) 
Jd 

which will be rederived below by a more direct method. 
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Corrections to the leading term. Taking into account the discrete "atomic" struc- 
ture of the Dyson gas, one is able to find the subleading corrections to the free energy. 

The first correction comes from a more accurate integral representation of the sum 
log \ zi — Zj\, when passing to the continuous theory. Namely, one should exclude 
the terms with i — j, writing 

lo S \zi-Zj\=^2 log \zi-z 3 \-J2 lo S \t( z i) I 

where t is a short-distance cutoff (which may depend on the point Zj). It is natural to 
take the cutoff to be 

«*> ~ \S (4 - 22) 

which is the mean distance between the charges around the point z. (In the context 
of the quantum Hall effect, £ ~ \fhjB is called the magnetic length.) This gives the 
improved estimate for E[p c i\. 

-{3h 2 E[ Pd ] = F + (3h J p cl (z)^p~ d Jz)d 2 z-^N\ogh + a l N (4.23) 

where oti is a numerical constant which can not be determined by this argument. 

Another correction comes from the integration measure when one passes from the 
integration over Zj to the integration over macroscopic densities 1 . We can write 

J] d 2 z 3 = N\ J[p] [Dp] 

j 

where [Dp] is an integration measure in the space of densities, J[p] is the Jacobian 
of this change of variables and the factor N\ takes into account the symmetry under 
permutations (all the states that differ by a permutation of the charges are identical). 
To estimate the Jacobian, we divide the plane into iV microscopic "cells" such that j-th 
particle occupies a cell of size £(zj), where £(zj) is the mean distance (4.22) between the 
particles around the point Zj. All the microscopic states in which the particles remain in 
their cells are macroscopically indistinguishable. Given a macroscopic density p, J[p] is 
then approximately equal to the integral f cejls Y[j d 2 Zj, with each particle being confined 
to its own cell. Therefore, J[p] ~ Ylj^ 2 ( z j)^ an d thus log J[p\ (sometimes referred to as 
entropy of the state with the macroscopic density p) is given by 

logJ[p] = ~\J Pd(z) log p d (z)d 2 z + Nlogh + a 2 N (4.24) 

where a 2 is a numerical constant. This result agrees with the corresponding Jacobian 
obtained within the collective field theory approach [34]. 

Combining (4.23), (4.24) with p = p d , and taking into account the factor iV! in the 
measure, we obtain: 

N 

c{N) = \ogN\ + —(2-(3)\ogh + aN (4.25) 



1 I thank A.Abanov for a discussion on this point. 
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where a is a numerical constant, and 



F 1/2 = - J a(z) log(na(z)) d 2 z (4.26) 

The term F 1 / 2 is thus the sum of the contribution due to the short- distance cutoff and 
the entropy contribution, which cancel each other in the ensemble of normal self-dual 
matrices (at (3 = 2). Another remarkable case when Fi/ 2 vanishes exactly 2 is the case of 
quasiharmonic potentials. 

The result for Fy 2 can be derived in a more rigorous way from the loop equation (see 
Appendix C). Here we simply note that the variation of (4.26) over the potential W does 
yield the correction to the mean density given by (4.16). One can verify this using the 
variational technique presented below in Section 4.4. 

The result for c(N) is in agreement with the dimensionality argument. It is easy to 
see that e F °^ h carries the dimension of [length]^ , and higher terms are dimensionless. 
The dimension of Zjq is given by (4.3). Therefore, must carry the residual dimension 
[length] N ( 2 ~P\ which agrees with (4.25). For quasiharmonic potentials, the constant a 
can be found explicitly: a = log a/27t 3 , as is readily seen from the Gaussian case. 

To summarize, the asymptotic expansion of the partition function as K — > has the 
form 

Z N = N\h^ N e aN ew(^ + ^f + F 1 + J2 ^ 2 F k/2 \ (4.27) 

V fc>3 / 

where Fq and F\/ 2 are given by (4.20) and (4.26) respectively. The higher corrections are 
due to fluctuations of the eigenvalues around the equilibrium configuration. No simple 
method to find their explicit form is known. In principle, these corrections can be found 
by expanding the loop equation in powers of h (see Appendix C), similarly to how it 
goes for the model of one Hermitian matrix [35, 36]. However, the calculations are rather 
tedious, even in the first two orders. At present only fragmentary results are available. 
Some of them look quite suggestive. 

For example, the result for the Fi-correction obtained in [37] for the case of the normal 
matrix model with quasiharmonic potential and a connected support of eigenvalues has 
a clean interpretation as the free energy of the theory of free bosons in the domain D c 
with the Dirichlet boundary conditions. Namely, the calculations yield the result 




(4.28) 



where z(w) is the univalent conformal map from D c onto the exterior of the unit cir- 
cle. Comparison with the Polyakov-Alvarez formula [38, 39] allows one to identify this 
quantity with — \ logdet(— A D c), where det(— A D c) is the regularized determinant of the 
Laplace operator in D c with the Dirichlet boundary conditions. This suggests the inter- 
pretation through free bosons 3 . Presumably, the higher corrections to the free energy are 
connected with spectral geometry of the Laplace operator, too. Recently, some progress 
in computation of F\ and, more generally, in understanding the structure of the whole 

2 However, its variational derivative, pn, does not vanish, as is seen from (4.16). 

3 A similar interpretation of the ^-correction in the Hermitian matrix ensemble was suggested in [40]. 
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Figure 4: Tie normal displacement of the boundary. 

series (including the case of disconnected supports) in models of Hermitian matrices was 
achieved [41]. Conjecturally, the answer is to be expressed in terms of a (conformal?) 
field theory on the complex curve T introduced at the end of Section 4.2. 

Finally, we note that the structure of the loop equation suggests to rearrange the 
/i-expansion of the free energy and to write it in the "topological" form F = J2 g >o ^ 9 F gi 

where each term has its own expansion in e = (2 — f3)h: F g = F® + Yl n >i ^Fg™' '• 



4.4 Correlation functions in the large N limit 

Variational technique and the Dirichlet boundary value problem. Correlation 
functions in the leading order in h can be obtained from the free energy by variation 
w.r.t. W(z) according to the formulas from Section 3.1. For a variation of the potential, 
W — > W + 5W with Nh fixed we ask how D changes. It is convenient to describe small 
deformations D — > D, by the normal displacement 6n(£) of the boundary at a boundary 
point f (Fig. 4). 

Consider a small variation of the potential W in the condition (4.10), which determines 
the shape of D at a fixed t. To take into account the deformation of the domain, 5D = 
D \ D, we write, for any fixed function /, 

s((f f(()d() = I /(c)dc = 2i/ df(()d 2 (^2t(f 8f(0H0\d(\ 

\JdD J Jd(SD) JSD JdD 

and thus obtain from (4.10): 

t amoK + j i AinoMQ m = o (4 . 29) 

JdD z ~ C 2 J dD z — Q 

Here the first term comes from the variation of W and the second one comes from the 
change of D. 

This is an integral equation for the Sn((). It can be solved in terms of the exterior 
Dirichlet boundary value problem. Given any smooth function f(z), let f H (z) be its 
harmonic continuation (already introduced in Section 4.2) from the boundary of D to its 
exterior, i.e., the function such that it is harmonic in D c , Af H = 0, and regular at oo, 
and f H (z) = f(z) for all z G dD. The harmonic continuation is known to be unique. 
Explicitly, a harmonic function can be reconstructed from its boundary value by means 
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of the Dirichlet formula 

f H {z) = ~i f(0d n G(z,0m (4.30) 

^ JdD 

The main ingredient of this formula is G(z,£), which is the Green function of the domain 
D c : 

A z G(z,() = 2ir5(z-() in D c , G(z, C) = if z G dD 

As C — » z, it has the logarithmic singularity G{z,C) — > log \z — (\. The well known 
properties of harmonic functions imply that d n G(z,£) < for all £ G <9D (the operator 
of the normal derivative acts here to the second argument). 

Consider the integral 

r d(5W H )d( 

JdD Z-C 

which is obviously equal to for all z inside D, subtract it from the first term in (4.29) 
and rewrite the latter in the form 

d(5W-5W H )d( _ i_ r d n (5W-5W H ) 

dD z ~ C 2 Jg D Z — ( 

Here the integral over dC, is transformed to the integral over the line element \d(\. The 
normal derivative is taken in the exterior of the boundary. After this simple transforma- 
tion our condition acquires the form 

^ l(C) + A(C) Kl = for all z e D (4.31) 

dD z ~ C 

where 8n\{z) := AW(z)Sn(z) (just for brevity) and R is the Neumann jump operator. 
Acting on a smooth function /, this operator gives the difference between the normal 
derivative of this function and the normal derivative of its harmonic extension: 

Rf(z)=d-(f(z)-f H (z)) (4.32) 

The superscript indicates that the derivative is taken in the exterior of the boundary. 

By properties of Cauchy integrals, it follows from (4.31) that Sni(z) + -R5W^(z)| 

is the boundary value of an analytic function h(z) in D c such that h(oo) = 0. For z G D c , 
this function is given by 

m = J_I fa l(C ) + fl(C) 

If h is not identically zero, the number of zeros of g in D c , counted with multiplicities, 
is given by the contour integral § dD d(argh). Since argh(z) = axg(\dz\/dz) = —0{z), 
where 9 is the angle between the tangent vector and the real axis, this integral is equal 
to 1. We conclude that the function h has exactly one simple zero outside the domain 
D. It is just the zero at oo. 

On the other hand, the variation of the normalization condition yields, in a similar 
manner: 

(Sri! + d n 5W)\d(\ = (4.33) 

dD 
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This relation implies that the zero at oo is at least of the 2-nd order. Indeed, expanding 
the Cauchy integral around oo, 

2mh(z) -+-<f (5n 1 +d n 5W-d n 5W H )\d(\ + 0( y z- 2 ) 

z JdD 

one concludes, using the Gauss law <f dD d n 5W H \dC,\ = — J DC ASW H d 2 ( = 0, that the 
coefficient in front of 1/z vanishes. We have got a contradiction. 

Therefore, h(z) = 0, and so Sn^z) + RSW(z) = 0. This gives the following result 
for the normal displacement of the boundary caused by small changes of the potential 
W -> W + 5W: 

* (w $;r w) (4 - 34) 



Some results for the correlation functions. In order to find the correlation func- 
tions of traces, we use the general variational formulas (3.5), where the exact free energy 
is replaced by the leading contribution (4.20): 

& {p{z)) = wU = pd{z) ' & m*)a»(*)>c = 

Basically, these are linear response relations used in the Coulomb gas theory [42]. In this 
approximation, the eigenvalue plasma is represented as a continuous charged fluid, so 
the information about its discrete microscopic structure is lost. So, these formulas give 
"smoothed" correlation functions in the first non-vanishing order in h. They are correct 
at distances much larger than the mean distance between the charges. 

Here are the main results for the correlation functions obtained by the variational 
technique. For details of the derivation see [7] and Appendix C. 

The leading contribution to the one-trace function was already found in Section 4.2. 
Here we present the general result including the first subleading correction which can be 
found by variation of (4.26): 

/?<tr/(S)> = ±J^a(z)f(z)d 2 z 

(4.35) 



8tt 



(l+\oga(z))Af(z)d 2 z- I \o Z a{z)Rj{z)\dz\ 

JdD 



+ 0(h) 



where R is the Neumann jump operator (4.32). Applying this formula to the function 
(p(z), we get the familiar result ((p(z)) = <fd( z ) + ¥n( z ) + 0{h 2 ), where fn, is given by 
(4.15). The connected two-trace function is: 

(3 (trf tig) c =^J d VfVgd 2 z - i- j fd n g H \dz\ + 0(h) (4.36) 

In particular, for the connected correlation functions of the fields (fi(zi), <p(z 2 ) (see (4.2)) 
this formula gives (if 2 12 G D c ): 

! (cp(z 1 )< f (z 2 )) c = G(z u z 2 ) - G(z u oo) - G(oo, z 2 ) - log ^ ~ ^ + 0(h) (4.37) 



2ph 
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where G is the Green function of the Dirichlet boundary value problem and 



cxp 



lim(log|f| + G(f,oo)) (4.38) 



is the (external) conformal radius of the domain D. The 2-trace functions are universal, 
i.e., they depend on the shape of the support of eigenvalues only and do not depend 
on the potential W explicitly. They resemble the two-point functions of the Hermitian 
2- matrix model found in [18]; they were also obtained in [43] in the study of thermal 
fluctuations of a confined 2D Coulomb gas. The structure of the formulas indicates that 
there are local correlations in the bulk as well as strong long range correlations at the 
edge of the support of eigenvalues. (See [44] for a similar result in the context of classical 
Coulomb systems). 

From the mathematical point of view, the significance of formula (4.37) is to provide a 
link between such seemingly unrelated disciplines as classical analysis in two dimensions 
and the random matrix theory. Namely, different limits or certain specifications of the 
arguments in this formula allow one to represent some important objects of classical 
analysis associated with the domains D and D c (e.g., the conformal map onto the unit 
circle and its Schwarzian derivative, the Bergman kernel) in terms of correlations between 
eigenvalues of random matrices. 

Further variation of the pair density correlation function suggests that, starting from 
n — 3, the connected n-point density correlations vanish in the bulk in all orders of H 
(in fact they are exponential in 1/h). The entire leading contribution comes from the 
boundary. The result for the connected three-trace function is: 

H f \dz 




16tt 2 



/ ^\f[RfAz) + 0(h 2 ) (4.39) 
JdD v{z) fJi 



5 The matrix model as a growth problem 
5.1 Growth of the support of eigenvalues 

When N increases at a fixed potential W, one may say that the support of eigenvalues 
grows. More precisely, we are going to find how the shape of the support of eigenvalues 
changes under t — > t + St, where t = Nh, if W stays fixed. 

The starting point is the same as for the variations of the potential, and the calcula- 
tions are very similar as well. Variation of the condition (4.10) and of the normalization 
condition yields 

r Aw(t)5n(t) = Q (for all z ^ D ^ r Aw ^ 5n ^ = _ 47r(36t 

JdD z ~ C JdD 

The first equation means that AW(z)5n(z)^j- is the boundary value of an analytic 
function h(z) such that h(z) = —4ir(35t/z + 0(z~ 2 ) as z — > oo. The solution for the 
6n(z) is: 

5n(z) = --^-d n G(oo,z) (5.1) 
Znayz ) 



33 



z(w) 



z-plane 



w-plane 



Figure 5: The mutually inverse conformal maps w(z) and z(w). 

where G is the Green function of the Dirichlet boundary problem in D c . For quasihar- 
monic potentials (with o — 1/ir), the formula simplifies: 

Sn(z) = ~5td n G(oo,z) (5.2) 

Identifying t with time, one can say that the normal velocity of the boundary V n = 5n/St, 
is proportional to gradient of the Green function: V n oc — d n G(oo, z). This result is quite 
general. It holds for any (not necessarily connected) domains of eigenvalues with a smooth 
boundary. 

If the domain is connected, the Green function can be expressed through the conformal 
map w(z) from D c onto the exterior of the unit circle: 



G(z 1 ,z 2 ) = log 



w(zi) - w(z 2 ) 



1 - w(z 1 )w(z 2 ) 



(5.3) 



In particular, G(oo,z) = — \og\w(z)\. As \z\ — > oo, w(z) = z/r + 0(1), where r is the 
external conformal radius of the domain D which enters eq. (4.37). It is easy to see that 
<9 n log \ w(z)\ = \w'(z)\ on <9D, so one can rewrite the growth law (5.2) as follows: 

5n(z) = ^5t\w'(z)\ (5.4) 

It is worth noting that the inverse map, z(w), is the classical (h — > 0) limit of the 
Lax operator (3.17) of the 2D Toda hierarchy. Indeed, making the rescaling nh = t, 
d/dn = hd t , we see that the shift operator e Wt can be replaced by a commuting variable 
w with the Poisson bracket {logu>, t} = 1. In this limit, known also as the dispersionless 
limit of the 2D Toda hierarchy, the Lax operator L(e Wt ) converts into the function z(w) 
given by the series of the form 

z(w) = r(t)w + u k(f)w~ k (5-5) 

k>0 

It defines the one-to-one conformal map from the exterior of the unit circle onto D c , 
which is inverse to the map w(z) (Fig. 5). For more details see [8]. 
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The basic formula which allows one to find the t-derivative of any quantity of the 
form J* D f(z) d 2 z (where the function / is assumed to be independent of t) immediately 
follows from (5.1): 

5 Qf f(z) d 2 z^j = £ f(z)5n(z) \dz\ =-^j M. dnG ^ z ) \ dz \ 

In the r.h.s. we recognize the value at infinity of the harmonic continuation of the function 
f(z)/a(z), so the result is 

§t(f D f(*)**) =P(f/°f (oo) (5.6) 

Using this formula and the integral representation of F , we find t-derivatives of the free 
energy. The first derivative, 

d t F = 2 a(z) log \z\ d 2 z (5.7) 

J D 

was already found in Section 4.3 by other means (see (4.21)). The second derivative is 
proportional to the logarithm of the conformal radius (4.38): 

= 2(3 lim (log |*| + G(z, oo)) = 2(3\ogr (5.8) 

Z = OQ Z^rOO 



d 2 t F, = 2P{\og\z\) 



In this connection let us also mention the nice formula for the conformal map w(z), 

w{z) = hm — — 5.9 

N^oo 1p N (z) 

which follows from the relation log (Pn+i(z)/Pn(z)) = Hd t (trlog(2; — $)) + 0(h) after 
calculating the t-derivative according to the above rule. This connection between con- 
formal maps and orthogonal polynomials goes back to the classical theory of analytic 
functions (see e.g. [45]). It is the context of the random matrix theory where it looks 
natural and easily understandable. 



5.2 Laplacian growth 

The growth law (5.2) is common to many important problems in physics. The class 
of growth processes, in which dynamics of a moving front (an interface) between two 
distinct phases is driven by a harmonic scalar field is known under the name Laplacian 
growth. The most known examples are viscous flows in the Hele-Shaw cell (the Saffman- 
Taylor problem), filtration processes in porous media, electrodeposition and solidification 
of undercooled liquids. A comprehensive list of relevant papers published prior to 1998 
can be found in [46]. Recently, the Laplacian growth mechanism was recognized [10] in 
a purely quantum evolution of semiclassical electronic blobs in the QH regime. 
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Figure 6: The Hele-Shaw cell. 

The Saffman- Taylor problem. Let us describe the main features of the Laplacian 
growth on the example of viscous flows. To be specific, we shall speak about an interface 
between two incompressible fluids with very different viscosities on the plane (say, oil 
and water). In practice, the 2D geometry is realized in the Hele-Shaw cell - a narrow 
gap between two parallel glass plates. In this version, the problem is also known as the 
Saffman- Taylor problem or viscous fingering. For a review, see [14]. The velocity field in 
a viscous fluid in the Hele-Shaw cell is proportional to the gradient of pressure p (Darcy's 
law): 

b 2 

V = -KV Pl K= — 

12fi 

Here the constant K is called the filtration coefficient, \x is viscosity and b is the size of 
the gap between the two plates. Note that if /i — > 0, then Vp — > 0, i.e., pressure in a 
fluid with negligibly small viscosity is uniform. Incompressibility of the fluids (VV = 0) 
implies that the pressure field is harmonic: Ap = 0. By continuity, the velocity of the 
interface between the two fluids is proportional to the normal derivative of the pressure 
field on the boundary: V n = —Kd n p. 

To be definite, we assume that the Hele-Shaw cell contains a bounded droplet of 
water surrounded by an infinite "sea" of oil (another possible experimental set-up is an 
air bubble surrounded by water). Water is injected into the droplet while oil is withdrawn 
at infinity at a constant rate, as is shown schematically in Fig. 6. The latter means that 
the pressure field behaves as p oc — log|z| at large distances. We also assume that the 
interface between oil and water is a smooth closed curve 7 which depends on time. As it 
was mentioned above, if viscosity of water is negligible, then one may set p = inside the 
water droplet. However, pressure usually has a jump across the interface, so p in general 
does not tend to zero if one approaches the boundary from outside. This effect is due 
to surface tension. It is hard to give realistic estimates of the surface tension effect from 
first principles, so one often employs certain ad hoc assumptions. The most popular one 
is to say that the pressure jump is proportional to the local curvature of the interface. 
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To summarize, the mathematical setting of the Saffman- Taylor problem is as follows: 



V n = -d n p on 7 

Ap = in oil 

p^ — log \z\ in oil as z — > oo (5.10) 

p = in water 

p(+) — p(-) — — VK across 7 

Here v is the surface tension coefficient and k is the local curvature of the interface. (The 
filtration coefficient is set to be 1.) The experimental evidence suggests that when surface 
tension is small enough, the dynamics becomes unstable. Any initial domain develops 
an unstable fingering pattern. The fingers split into new ones, and after a long lapse of 
time the water droplet attains a fractal-like structure. This phenomenon is similar to the 
formation of fractal patterns in the diffusion-limited aggregation. 

Comparing (5.2) and (5.10), we identify the D and D c with the domains occupied 
by water and oil respectively, and conclude that the growth laws are identical, with the 
pressure field being given by the Green function: p(z) = G(oo,z), and p = on the 
interface. The latter means that supports of eigenvalues grow according to (5.10) with 
zero surface tension, i.e., with v = in (5.10). 

Neglecting the surface tension effects, one obtains a good approximation unless the 
curvature of the interface becomes large. We see that the idealized Laplacian growth 
problem, i.e., the one with zero surface tension, is mathematically equivalent to the 
growth of the support of eigenvalues in ensembles of random matrices N, !N° and C. 
This fact clarifies the origin of the integrable structure of the Laplacian growth with 
zero surface tension discovered in [47]. The link to the normal matrix model has been 
established in [6], see also [48]. 



The finite-time singularities. As a matter of fact, the Laplacian growth problem 
with zero surface tension is ill-posed since an initially smooth interface often becomes 
singular in the process of evolution, and the solution blows up. The role of surface tension 
is to inhibit a limitless increase of the interface curvature. In the absense of such a cutoff, 
the tip of the most rapidly growing finger typically grows to a singularity (a cusp). In 
particular, a singularity necessarily occurs for any initial interface that is the image of 
the unit circle under a rational conformal map, with the only exception of an ellipse. 

An important fact is that the cusp-like singularity occurs at a finite time t = t c , i.e., 
at a finite area of the droplet. It can be shown that the conformal radius of the droplet 
r (as well as some other geometric parameters), as t — > t c , exhibits a singular behaviour 

r - r c oc (t c - ty 1 

characterized by a critical exponent 7. The generic singularity is the cusp (2,3), which 
in suitable local coordinates looks like y 2 = x 3 . In this case 7 = — |. The evolution can 
not be extended beyond t c . 

A similar phenomenon was well-known in the theory of random matrices for quite a 
long time, and in fact it was the key to their applications to 2D quantum gravity and 
string theory. In the large N limit, the random matrix models have critical points - the 
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points where the free energy is not analytic as a function of a coupling constant. As we 
have seen, the Laplacian growth time t should be identified with a coupling constant of 
the normal or complex matrix model. In a vicinity of a critical point, 

F ~F rcs +a (t c -t) 2 ^ 

where the critical index 7 (often denoted by 7 str in applications to string theory) depends 
on the type of the critical point. Accordingly, the singularities show up in correlation 
functions. Using the equivalence established above, we can say that the finite-time blow- 
up (a cusp-like singularity) of the Laplacian growth with zero surface tension is a critical 
point of the normal and complex matrix models. 

5.3 The semiclassical limit for electrons in magnetic field 

As we have seen in Section 2.3, the system of iV electrons in the plane in non-uniform 
magnetic field, which fully occupy the lowest energy level, is equivalent to the ensemble 
of normal N x N matrices, where N is degeneracy of the level. In this section we study 
the QH droplet in the semiclassical regime. The equivalence with the large N limit of 
the matrix model suggests to identify the semiclassical QH droplet with the support of 
eigenvalues. Remarkably, it is this limit where one makes contact with the purely classical 
Saff man- Taylor problem. 

A remark is in order. The limit S^Owe are talking about is really a semiclassical 
limit, or better to say "partially classical". Although the Planck constant h tends to 
zero, all the particles remain at the lowest (most quantum) energy level, assuming their 
mass is small or the magnetic field is large. The true quasiclassical limit would imply 
that the particles occupy higher energy levels. 

We recall that the joint probability to find electrons at the points Zi is \^n(zi, ■ ■ ■ , zn)\ 2 , 
with oc detjVxAr|/0n(zfc)]' where 

iMz) = -^=P n - 1 {z)e w ^/^ (5.11) 

are orthogonal one-particle wave functions for electrons in the magnetic field B = — |AW 
at the lowest level E = 0. The level is assumed to be completely filled, i.e., n = 
0,1 . . . , N = [<f)/4>o\, where 0o is the flux quantum. Then the mean density of the elec- 
trons coincides with the expectation value of the density of eigenvalues in the normal or 
complex matrix model. 

The degeneracy of the level cam be controlled in different ways. One of them is 
to assume the following arrangement (see Fig. 7). Let a strong uniform magnetic field 
Bq > be applied in a large disk of radius Rq. The disk is surrounded by a large annulus 
Rq <\z\ < R\ with a magnetic field B\ < such that the total magnetic flux through the 
system is N(p : tiBqR^— nBi(R\ — R^) = N(p Q . The magnetic field outside the largest disk 
\z\ < R\ vanishes. The disk is connected through a tunnel barrier to a large capacitor 
that maintains a small positive chemical potential slightly above the zero energy. If the 
field Bq is strong enough, the gap between the energy levels is large, and the higher levels 
can be neglected. In the case of the uniform fields B and B\ the QH droplet is a disk of 
radius r = \J Nh <C R trapped at the origin. 
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the droplet 



Figure 7: The configuration of magnetic fields. 

In this set-up, let us apply a non-uniform magnetic field, SB, somewhere inside the 
disk \z\ < Ro but well away from the droplet. This leads to the following two effects. 

The Aharonov-Bohm effect in the QH regime [10, 49]. Suppose that the nonuni- 
form magnetic field SB does not change the total flux: J 5Bd 2 z = 0. As is argued above, 
the shape of the droplet is the same as that of the support of eigenvalues in the ensemble 
N with the potential 

w(z) = -^\A 2 - l -J log \z - ci mo d 2 c 

The second term is harmonic inside and around the droplet. One may have in mind 
thin solenoids carrying magnetic flux ("magnetic impurities"). In the case of point-like 
magnetic fluxes at points Oj the change of the potential is SW(z) = qi log \ z — Oj|. 

Let us stress that in the presence of the fluxes, the shape of the droplet is no longer 
circular although the magnetic field inside the droplet and not far from it remains uniform 
and is not changed at all (Fig. 8). In this respect this phenomenon is similar to the 
Aharonov-Bohm effect. Due to the quantum interference the electronic fluid is attracted 
to positive fluxes and is repelled by negative ones. The response of the droplet to an 
infinitesimal change of the magnetic field SB is described by eq. (4.34) in which 

5W H {z) - SW(z) = - [ G(z,()SB(()d 2 ( 

In fact this formula holds for arbitrary SB, not necessarily vanishing inside the droplet. In 
particular, for small point-like fluxes 5qi at some points a$ we have 5W = ^ Sqi log \z— a^j , 
SB = —nY^i$<lib { ' 2 \ z ~ a i)> an d SW H (z) — SW(z) = G(z,ai)Sqi. If a { is inside, 

G(z, di) is set to be zero. The sum, therefore, is over outside fluxes only. The fluxes inside 
the droplet, if any, appear to be completely screened and do not have any influence on 
its shape. 
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Figure 8: The electronic droplet in the presence of magnetic impurities. 

Growth of the electronic droplet. If the total magnetic flux increases, with magnetic 
impurities kept fixed, the electronic droplet grows. For example, one may adiabatically 
increase B\, with Bq and 5B fixed. Then the droplet grows because the degeneracy of the 
lowest level is enlarged and new electrons enter the system. The growth is described by 
eq. (5.1) with AW(z) = —2B which is equivalent to the Darcy law. This phenomenon is 
purely quantum. Like the Aharonov-Bohm effect, it is caused by quantum interference. 
Its characteristic scale is less than that of the Saff man- Taylor fingering by a factor of 
10 9 . The correspondence established above suggests that the edge of the QH droplet 
may develop unstable features similar to the fingers in the Hele-Shaw cell. 

5.4 Semiclassical ^-function 

In this Subsection we derive the semiclassical asymptotics of the ^-function (5.11), i.e., 

1pN+l{z) ^ l/>(z) 

To avoid cumbersome technical details, we mainly consider models with quasiharmonic 
potential. 

First of all, it is necessary to know the large N limit of the orthogonal polynomials 
P/v- Truncating the general formula (3.6) at the second term in the exponent, we can 
write: 

P N ( Z ) = (det(z - $)) = (J e tr iog(z-*)^ 

(5.12) 

h = exp (V log(z - $)> + \ ((tr log(* - $)) 2 > c + . . 

The first term in the r.h.s. is 0(h~ 1 ), the second one is 0(h°) and the ignored terms 
vanish as h — > 0. However, this formula should be applied with some care since the 
logarithm is not a single-valued function. This formula is correct if one can fix a single- 
valued branch of the logarithm. It is possible if z is outside D. For z inside D an analytic 
continuation should be used. To take care of the normalization, we also need 

h N = {n Z "^ Zn = {^h) l/2 exp (h-^Fo + \d 2 t F Q + . . .) (5.13) 

which is written here with the same precision as (5.12). (We have used (4.27) with 
a = log V27Y 3 h and have taken into account that Fi/ 2 = for quasiharmonic potentials.) 
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Let us first keep the dominant terms in the r.h.s. of (5.12), (5.13) and ignore the 
O(h ) terms for a while. We have, for z E D c : \4>(z)\ 2 ~ e ^W-<vW» > r, using 

the results of the previous section, 

\ij(z)\ 2 ~ e- 2 " 4 ^, = ^(z) - ¥>d(0) - 

As we know, A(z) defined by this formula is zero on the boundary. The analytic contin- 
uation of the A(z) inside the contour 7 = <9D can be done using the Schwarz function. 
Namely, since 2dA(z) = z — S(z), the desired analytic continuation can be defined (up 
to a constant) by the formula 



2A(z) = \z\ 2 - 27Ze j S(() d( 

Clearly, A(z) defined in this way is constant on the contour 7. Indeed, if z ij2 € 7, then 

2(A(z 2 ) - A(zi)) = \z 2 \ 2 - \z x \ 2 - 2TZe I S(z)dz 

J 21 



' 21 

pZ2 PZ2 rZ2 

= z 2 z 2 — Z1Z1 — / zdz — / zdz = z 2 z 2 — z{z\ — \ d(zz) = 

J Zl J Zl J Zl 

Since A(z) should be zero on 7, we finally define 

A(z) = \\z\ 2 - \U 2 -Kej*S(()d( (5.14) 

where £0 is an arbitrary point on the contour. We call the function defined by (5.14) 
the effective action. From the above it follows that its first derivatives vanish for all 
z £ 7: dA(z) = dA(z) = 0. This means that \4>\ 2 ~ e^ 2 ^^ has a sharp maximum on the 
contour 7. We may say that purely quantum particles are in general delocalized in the 
plane, purely classical particles are localized at some points in the plane while partially 
classical particles, like our electrons, are localized on closed curves in the plane. 

Let us turn to the O(ft )-corrections, which give the subexponential factor in the 
asymptotics of the ^-function 4 . We assume, for simplicity, that the domain D is con- 
nected. Extracting analytic and anti- analytic parts of eq. (4.37), we get 

<(tr log(z - $)) 2 ) c = \og{rw\z)) + 0(h) 

where r is the external conformal radius of D and w'(z) is the derivative of the conformal 
map from D c onto the exterior of the unit circle. Plugging this into (5.12) and taking 
into account eq. (5.8), we finally obtain: 

l^)| 2 = ^=e- 2 ^)/* (5.15) 

This formula does resemble the WKB asymptotics in quantum mechanics. If a singularity 
of the conformal map is sufficiently close to the boundary from inside, the asymptotics 
becomes invalid in this region. 

4 In what follows we restrict our consideration to the squared modulus ^(z)! 2 . The -0-function itself 
contains a phase factor which rapidly oscillates in both tangential and normal directions. (I thank 
A.Boyarsky and O.Ruchaisky for a discussion on this point.) 
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The effective action can be expanded near the contour, where it takes the minimal 
value: 

A(z + 5 n z) = \5 n z\ 2 t -K(z)\5 n z\ 3 + -^(z)^^ + ... (5.16) 

Here k(z) is the local curvature of the contour at the point z and 5 n z is a small deviation 
from the point z G 7 in the normal direction. (The upper and lower signs correspond to 
the outward and inward deviations respectively.) A similar expansion of log re& ds 

log \w'(z + 5 n z)\ = log \w'(z)\ ± (\w'(z)\ - k{z)) \8 n z\ + ... (5.17) 

Some details of the derivation are given in Appendix E. Therefore, if n(z) <C ft _1//2 , the 
squared modulus of the ^-function is well approximated by the sharp Gaussian distribu- 
tion in the normal direction with the amplitude slowly modulated along the curve: 

\^(z + 5 n z)\ 2 ~Pe" 2 M 2 / s (5.18) 
V27T 3 /j 

In the case of general potentials the calculations lead to a similar result: 



\ij{z + 5 n z)\ 2 ~ J^±\ w \ z )\e- 2 ^\ & - z \ 2 ' h (5.19) 

We see that the width of the Gaussian distribution depends on the point of the curve 
through the function a (which is proportional to the magnetic field in the QH interpreta- 
tion). Note that this asymptotics is consistent with the normalization J \ifj(z)\ 2 d 2 z = 1. 
The easiest way to see this is to notice that the Gaussian function in (5.18) (as well as 
the one in (5.19)) tends to the delta function 8{z;^) with the support on the curve (see 
Appendix B). Hence one can formally write the limiting ^-function in the form 

mz)\ 2 = l -^5(z ]1 ) (5.20) 

and the right normalization is transparent. 

Finally, we note that the growth law (5.4) (the Darcy law) can be written in the 
suggestive form 

V n (z)^\^(z)\ 2 (5.21) 
The normal velocity V n is defined by the relation 

d t Q(z;D(t)) = VJ(z; dD(t)) 

where Q(z; D) is the characteristic function of the domain D. Since the classical value of 
(p(z)) is Q(z; D)/n, we can represent the Darcy law in yet another form: 

o t (p(z)) = mz)\ 2 

In fact it is the h — > limit of the exact relation 

(p(z)) N -(p(z)} N _ 1 = h\^ N (z)\ 2 (5.22) 

which immediately follows from the fact that {p{z)) N = hK^(z, z) and from the definition 
of the kernel function (3.16). In fact it is equivalent to the Hirota equation (3.12). We 
see that (5.22) can be regarded as a "quantization" of the Darcy law. 
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Appendices 
Appendix A 



The proof of the determinant representation (3.7). It is simpler to start with 
detCy, where is given in (3.8), and to show that AHdet coincides with Z N . The 
proof is a chain of obvious equalities. By definition, 

NldetCij = ^(-) Q (-) P Cq(i) > p ( i) . . . C Q ( N ),p(n) 
Q,p 

Plugging the explicit form of the CV,, we write: 

r N 

NldetCij = E(-) Q (-) P /n [*? W ~ 1 < (0 ~ 1 e WW d 2 ; 

Q,P J i=l 

and, interchanging the order of summation and integration, obtain the result: 



N\ det dj = 



N 



b-> q it.' 

. Q 



Q(i)-1 



i=l 



N 



p(o-i 



i=l 



N 



l[e w ™<Pz k = Z N 



k=i 



A]v(2i) 



Orthogonality of the polynomials (3.13). The definition (3.13) implies 

1 r n 



Z 3 



i=i 

We want to show that these polynomials are orthogonal in the complex plane. 

Set e w ^d 2 z = d\x for brevity. It is enough to show that J P n (z)z m dfi = for all 
m < n, i.e., 

f dfi(z)z r 

3 = 

Note that A n (zi, ...,z n ) Yl" =1 (z - Zj) = A n+1 (z 1 , ...,z n ,z). Setting z = z n+1 , we have: 

n+1 



/n 
\A n (z i )\ 2 ]l(z-z j )dfi(z j ) = 

7=1 



/it~r± 
A n+1 (^)A^)C+i Udfiizj) 

7=1 
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j n+1 » n+1 

= 7— Y^(-l)' +n+1 / A n+1 (z!, . . . , Z n+1 ) A n ( Zl , . . . , J h . . . , Zn +1 )z? H dn( Zj ) 
U 1=1 J j=l 

One can notice that the summation gives the expansion of the determinant 

1 ~ ~n—l ~m 

1 Z\ . . . Z^ Z-y 

1 ~ ~n—l ~m 

1 Z2 ■ ■ ■ Z 2 Z 2 

1 ~ ~n— 1 m 

1 ^ n+1 . . . Z n+1 Z n+1 

under the bar. If m < n, it vanishes. If m — n, it equals A n+1 (zi), and we get LHS = 

^l Z n+l- 

Appendix B 

Here we list some standard formulas often used in Section 4. 



The complex notation. 

• Complex coordinates: z = x + iy,z = x — iy,d z = \{d x — id y ), d z = \{d x + id y ). 

• The Laplace operator: A = d\ + dy = 4d z d s 

• Contour integrals: let / and g be any smooth functions defined in some neighbor- 
hood of the contour 7, then 



j> gd n f\dz\ = -2i & gdjdz -i j> fdg 



Singular functions. 



• Two-dimensional 5-function: S(z) = ^AlogN = -dz(l/z). The characteristic 
property of the delta-function is f(z)S(z — a) d 2 z = f (a) for any (smooth) function 
/■ 

• The 5-function with the support on a curve (a closed contour) 7: a function S(z; 7) 
such that 



J f(z)5(z; 1 )d 2 z = j f{z)\d 
for any smooth function /. 



• The "normal derivative" of the ^-function of the contour 7: a function 5'(z; 7) such 
that 



J f(z)5'(z; 1 )d 2 z = - jd n f{z)\dz\ 



"1 

for any smooth function /, with the normal vector being directed to the exterior 
of the contour. 

The characteristic function of the domain D: Q(z; D) = 1 if z G D and otherwise; 
VO(z;D) = -HS(z;dD). 
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Integral formulas. 

• Cauchy's integral formula (/ is any smooth function): 



1 I /(CR 1 f d- c f(()dX _ / -f( z ) , z E D 



2%i J dD z - 
In particular, f(()d( = 2% J D df(() d 2 (. 
• The Green formula: 



fAgd 2 z = - / VfVgd 2 z+ (t fd n g\dz\ 

D JD JdD 



where the normal vector looks outward D. 
The Dirichlet formula: 



u{z) = -^j>u{Qd n G{z,C)\dC\ 



for any function u harmonic in D c = C \ D. Here G(z, () is the Green function of 
the Dirichlet boundary value problem in D c . 



Appendix C 

Let us present some details of the ^-expansion of the loop equation (4.4). First of all we 
rewrite it in the form 



1 

2^ 



L(z,C) (A V (())d 2 ( = {dMz)Y - ((d( V (z)- Vd (z)) 2 ) - (2-0)h(^)) 



which is ready for the ^-expansion. Here 

dW{Q-dip cl {z) 



is the kernel of the integral operator in the l.h.s. (the "loop operator"). The zeroth order 
in h gives equation (4.7) which implies the familiar result <p c i(z) = — f D log \z — (\ 2 a(()d 2 ( 
for the ipd. To proceed, one should insert the series 

(ifi(z)) = <p d (z) + hip 1/2 (z) + h 2 Vl {z) + 0{h 3 ) 

(which corresponds to the /i-expansion (4.17) of the free energy) into the loop equation 
and separate terms of order h, h 2 etc. (In the notation adopted in the main body of the 
paper h(pi/ 2 — <fin + 0(h 2 ).) The terms of order h and h 2 give: 

± J L(z,C) (A^ 1/2 (())d 2 ( = -(2-P)d 2 Mz) 

L(z,0 (A Vl (0)d 2 C = - \{d Vl/2 (z)) 2 + (2-/?)0V 1/2 (z)l -u>(z) 



1 

2^ 
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where 



u(z) 



lim 



h~ 2 lim {d<p(z)<p(z')) < 



is the connected part of the pair correlator at merging points. If the point z is in D c , 
then eq. (4.37) yields 

(dtp(z) ip(z')) c = 2(3h 2 d z d z , (G(z, z') - log \z -z'\)+ 0(h 3 )) 

Since the r.h.s. is regular for all z, z' G D c , the points can be merged without any 
regularization and the result does not depend on the particular limit z' — > z. We thus 
obtain that the function oj{z) is proportional to the Schwarzian derivative of the conformal 
map w(z): 

pfw"'(z) 3 fw"{zy 2 " 
6 



u(z) 



w'(z) 



2 V w'(z 




The expansion of the loop equation can be continued order by order. In principle, this 
gives a recurrence procedure to determine the coefficients <Pk(z). However, the equations 
of the chain are integral equations in the plane, and it is not easy to solve them explicitly. 
Another difficulty is that in general one can not extend these equations to the interior of 
the support of eigenvalues because the ^-expansion may break down or change its form 
there. Indeed, in the domain filled by the gas of eigenvalues the microscopic structure of 
the gas becomes essential, and one needs to know correlation functions at small scales. 
Nevertheless, at least in the first two orders in h 2 the equations above can be solved 
assuming that z G D c . Note that in this region all the functions ipk(z) are harmonic. If 
these functions are known, the corresponding expansion coefficients of the free energy in 
(4.17) can be obtained by "integration" of the variational formulas (3.5). 

The procedure of solving the loop equation in the first two orders in h is too technical 
to be presented here. In the order h one is able to find a complete solution which gives 
formulas (4.15) and (4.26) mentioned in the main text. The solution in the next order, 
h 2 , is much more difficult to obtain. The results for </?i and F\ are still not available 
in full generality (i.e., for general (3 and W). Nevertheless, for normal matrices with a 
general potential {(5 = 1) and with a connected support of eigenvalues the iq-correction 
to the free energy can be found explicitly by the method outlined above. Here we present 
the result (mostly for an illustrative purpose), using the notation introduced in the main 
text (see (4.14), (5.5)): 



Fi = 



— |— / flog \z'(w)\d n log \z'(w)\ + 2 log \z\w)\)\dw\ 
J\ w \ =1 \ J 



+ 



24tt 
1 

8^ 



f \V X \ 2 d 2 z + 2(f K X \dz\ 

D JdD 

[ \V X \ 2 d 2 z- <f X dnX H \dz 

JD JdD 



1 



167T J D 



Axd 2 z + c 



where Co is a numerical constant and k(z) = d n log 
boundary. For quasiharmonic potentials only the first integral survives. 



w(z) 



w'(z) 



is the local curvature of the 
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Appendix D 

Here we demonstrate how the variational technique works with the 2-trace correlator 
(4.36). We shall use the variational formulas (3.5) in the following equivalent version. 
Set 5W(z) = eg(z), where g is an arbitrary smooth function and e — > 0. Then, in the 
first order in e, 

h5(trf) = e(trftrg) c 

This relation allows one to find the connected part of the two-trace correlation function 
by variation of the known one-trace function. Similar formulas hold for variations of 
multi-trace functions. 

We have, in the leading order in h: 

h(35 (tr /) = 5 [ I afd 2 z) = [ 5afd 2 z + [ afd 2 z = h + I 2 
\Jd J Jd Jsd 

The first integral, ii, can be transformed using the definition of a (4.9) and the Green 
formula: 

h = -±J B A(*K0/ * * = £ I V, V/ * z - ± £ fd„g W 
The second integral is 

h = -^f^W{z)f{z)8n{z)\dz\ 

where 5n(z) is to be taken from eq. (4.34): 5n(z) = ed n (g H (z)— g(z))/AW(z). Summing 
the two contributions, we get (4.36). 

Appendix E 

In this Appendix we obtain the expansions of the effective action A(z) 

Mz) = l -\z\ 2 - l -U 2 -Tie S(QdC 

near the contour 7 = <9D. We know that the first variation of A(z) vanishes on 7. To 
find the second variation, we write 

2SA{z) = \Sz\ 2 -TZe{S'(z){Sz) 2 ), z e 7 

Let us represent 5z as a sum of normal and tangential deviations w.r.t. the curve: 
5z = 5 n z + 8 t z, then 

2-Re {S'{z){5z) 2 ) = Tie {S'{z){5 n z) 2 ) + Tie {S'{z){5 t z) 2 ) + Tie {S'{z)5 n z5 t z) 

Using the obvious relations 

S t z IJ^z __ 1 
|^| ~ V S t z ~ y/S>(zj 



S n z 



5 t z 



5 t z 
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where the upper (lower) sign should be taken for the outward (inward) deviation, and 
the formula for the scalar product of 2-vectors x, y represented as complex numbers x, y, 
(x, y) = IZe (xy), we see that the first and second terms in the r.h.s. are equal to — \S n z\ 2 
and \S t z\ 2 respectively while the third one vanishes since the vectors 5 n z and 5 t z are 
orthogonal. Since \Sz\ 2 = \S n z\ 2 + |<5tz| 2 , we obtain the desired result 5A(z) = \5 n z\ 2 . 

The next terms of the expansion of the A(z) around the contour can be found in a 
similar way. They are expressed through the curvature k and its derivatives w.r.t. the 
arc length s along the curve. To perform the calculations in next two orders we need the 
following formulas for the k and k' = dn/ds through the Schwarz function [31]: 

( , * S"(z) „ , i S"'(z)S'(z) - l(S"(z)) 2 

K[Z) ~ 2 (S>(z))V 2 ' K [Z) ~ 2 (S'(z)f ' Z fc 7 

For z on the contour we have 

A(z + 5 n z) = \5 n z\ 2 - l -TZe (S"(z)(5 n z) 3 ) - ^JZe (S"\z)(5 n z) 4 ) + ... 
Now, with the help of the formulas for the curvature, it is easy to find that 

3/2 



4 



Tie (S"(z)(5 n z) 3 ) = 2K(z)7Ze (j (v^i)^^)^ = 2K(z)Ue ^ (^j (5 n 
whence 

Ke(S"(z)(5 n zf) = ±2K(z)\5 n z\ 3 

A similar computation gives IZe (S"'(z)(5 n z) A ) = — 6k 2 (z)\S ji z\ 4: and we obtain the expan- 
sion (5.16). The expansion of |u/(;z)| around the contour can be easily performed with 
the help of the relations 

w(z) 



k(z) = d n log 
valid for z G 7. 



w'(z) 



d n \og\w(z) \ = \w'(z) 
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